//
// This file is part of the Marble Virtual Globe.
//
// This program is free software licensed under the GNU LGPL. You can
// find a copy of this license in LICENSE.txt in the top directory of
// the source code.
//
// Copyright 2014 Gerhard Holtkamp
//

/* =========================================================================
  Procedures for calculating the positions of planets.
  The procedures in this unit are taken from Montenbruck, Pfleger
  "Astronomie mit dem Personal Computer", Springer Verlag, 1989
  and modified correspondingly.
  The library Astrolib has to be included.

  Copyright :Gerhard HOLTKAMP          26-MAR-2014   
  ========================================================================= */

#include "astr2lib.h"
#include <cmath>
using namespace std;

#include "astrolib.h"

extern double frac (double f);
extern double atan21 (double y, double x);


/*-------------------- Class Plan200 --------------------------------------*/
/*
 Ecliptic coordinates (in A.U.) and velocity (in A.U./day)
 of the Planets for Equinox of Date given in Julian centuries since J2000.
 ======================================================================
*/
Plan200::Plan200 ()
  { }

void Plan200::addthe (double c1, double s1, double c2, double s2, 							double& cc, double& ss)
{
 cc=c1*c2-s1*s2;
 ss=s1*c2+c1*s2;
}

void Plan200::term (int i1, int i, int it, double dlc, double dls, double drc,
			  double drs, double dbc, double dbs)
{
 if (it == 0) addthe (c3[i1],s3[i1],c[i],s[i],u,v);
 else
  {
   u=u*tt;
   v=v*tt;
  }
 dl = dl + dlc*u + dls*v;
 dr = dr + drc*u + drs*v;
 db = db + dbc*u + dbs*v;
}

Vec3 Plan200::velocity()   // return last calculated planet velocity
{
 return vp;
}

void Plan200::state (Vec3& rs, Vec3& vs)
{
 /* State vector rs (position) and vs (velocity) of the Sun in
    ecliptic of date coordinates for last calculated planet
 */
 rs = rp;
 vs = vp;
}

void Plan200::posvel ()
{
 /* auxiliary program to calculate position and velocity
    vectors rp, vp.
    to be called by the various planet procedures.
 */

 double	cl, sl, cb, sb;

 cl = cos(l); sl = sin(l); cb = cos(b); sb = sin(b);
 rp[0] = r*cl*cb; rp[1] = r*sl*cb; rp[2] = r*sb;

 // velocity
 vp[0] = (dr*cl*cb - dl*r*sl*cb - db*r*cl*sb) * 1E-4;
 vp[1] = (dr*sl*cb + dl*r*cl*cb - db*r*sl*sb) * 1E-4;
 vp[2] = (dr*sb                 + db*r*cb) * 1E-4;
}

/*--------------------- Mercury ------------------------------------------*/

Vec3 Plan200::Mercury (double t)   // position of Mercury at time t
{
 /*
   Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
   of Mercury for Equinox of Date.
   t is the time in Julian centuries since J2000.
 */

 const double arc = 206264.81;    // 3600*180/pi = ''/r
 const double p2 = M_PI * 2.0;

 int i;

 tt = t;
 dl = 0.0; dr = 0.0; db = 0.0;
 m1 = p2 * frac(0.4855407 + 415.2014314*t);
 m2 = p2 * frac(0.1394222 + 162.5490444*t);
 m3 = p2 * frac(0.9937861 + 99.9978139*t);
 m5 = p2 * frac(0.0558417 + 8.4298417*t);
 m6 = p2 * frac(0.8823333 + 3.3943333*t);
 c3[1] = 1.0; s3[1]= 0.0;
 c3[2] = cos(m1); s3[2] = sin(m1); c3[0] = c3[2]; s3[0] = -s3[2];
 for (i=3; i<11; i++)
	addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);

 // perturbations by Venus
 c[5] = 1.0; s[5] = 0.0; c[4] = cos(m2); s[4] = -sin(m2);
 for (i=4; i>0; i--)
	addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
 term(2, 5, 0, 259.74, 84547.39, -78342.34,0.01,11683.22,21203.79);
 term(2, 5, 1, 2.3, 5.04, -7.52, 0.02, 138.55, -71.01);
 term(2, 5, 2, 0.01, -0.01, 0.01, 0.01, -0.19, -0.54);
 term(3, 5, 0, -549.71, 10394.44, -7955.45, 0.0, 2390.29, 4306.79);
 term(3, 5, 1, -4.77, 8.97, -1.53, 0.0, 28.49, -14.18);
 term(3, 5, 2, 0.0, 0.0, 0.0, 0.0, -0.04, -0.11);
 term(4, 5, 0, -234.04, 1748.74, -1212.86, 0.0, 535.41, 984.33);
 term(4, 5, 1, -2.03, 3.48, -0.35, 0.0, 6.56, -2.91);
 term(5, 5, 0, -77.64, 332.63, -219.23, 0.0, 124.4, 237.03);
 term(5, 5, 1, -0.7, 1.1, -0.08, 0.0, 1.59, -0.59);
 term(6, 5, 0, -23.59, 67.28, -43.54, 0.0, 29.44, 58.77);
 term(6, 5, 1, -0.23, 0.32, -0.02, 0.0, 0.39, -0.11);
 term(7, 5, 0, -6.86, 14.06, -9.18, 0.0, 7.03, 14.84);
 term(7, 5, 1, -0.07, 0.09, -0.01, 0.0, 0.1, -0.02);
 term(8, 5, 0, -1.94, 2.98, -2.02, 0.0, 1.69, 3.8);
 term(9, 5, 0, -0.54, 0.63, -0.46, 0.0, 0.41, 0.98);
 term(10, 5, 0, -0.15, 0.13, -0.11, 0.0, 0.1, 0.25);
 term(0, 3, 0, -0.17, -0.06, -0.05, 0.14, -0.06, -0.07);
 term(1, 4, 0, 0.24, -0.16, -0.11, -0.16, 0.04, -0.01);
 term(1, 3, 0, -0.68, -0.25, -0.26, 0.73, -0.16, -0.18);
 term(1, 0, 0, 0.37, 0.08, 0.06, -0.28, 0.13, 0.12);
 term(2, 4, 0, 0.58, -0.41, 0.26, 0.36, 0.01, -0.01);
 term(2, 3, 0, -3.51, -1.23, 0.23, -0.63, -0.05, -0.06);
 term(2, 2, 0, 0.08, 0.53, -0.11, 0.04, 0.02, -0.09);
 term(2, 0, 0, 1.44, 0.31, 0.3, -1.39, 0.34, 0.29);
 term(3, 4, 0, 0.15, -0.11, 0.09, 0.12, 0.02, -0.04);
 term(3, 3, 0, -1.99, -0.68, 0.65, -1.91, -0.2, 0.03);
 term(3, 2, 0, -0.34, -1.28, 0.97, -0.26, 0.03, 0.03);
 term(3, 1, 0, -0.33, 0.35, -0.13, -0.13, -0.01, 0.0);
 term(3, 0, 0, 7.19, 1.56, -0.05, 0.12, 0.06, 0.05);
 term(4, 3, 0, -0.52, -0.18, 0.13, -0.39, -0.16, 0.03);
 term(4, 2, 0, -0.11, -0.42, 0.36, -0.1, -0.05, -0.05);
 term(4, 1, 0, -0.19, 0.22, -0.23, -0.2, -0.01, 0.02);
 term(4, 0, 0, 2.77, 0.49, -0.45, 2.56, 0.4, -0.12);
 term(5, 0, 0, 0.67, 0.12, -0.09, 0.47, 0.24, -0.08);
 term(6, 0, 0, 0.18, 0.03, -0.02, 0.12, 0.09, -0.03);

 // perturbations by Earth
 c[4] = cos(m3); s[4] = -sin(m3);
 for (i=4; i>1; i--)
	addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
 term(1, 1, 0, -0.11, -0.07, -0.08, 0.11, -0.02, -0.04);
 term(2, 4, 0, 0.1, -0.2, 0.15, 0.07, 0.0, 0.0);
 term(2, 3, 0, -0.35, 0.28, -0.13, -0.17, -0.01, 0.0);
 term(2, 1, 0, -0.67, -0.45, 0.0, 0.01, -0.01, -0.01);
 term(3, 3, 0, -0.2, 0.16, -0.16, -0.2, -0.01, 0.02);
 term(3, 2, 0, 0.13, -0.02, 0.02, 0.14, 0.01, 0.0);
 term(3, 1, 0, -0.33, -0.18, 0.17, -0.31, -0.04, 0.0);

 // perturbations by Jupiter
 c[4] = cos(m5); s[4] = -sin(m5);
 for (i=4; i>2; i--)
	addthe(c[i],s[i],c[4],s[4],c[i-1],s[i-1]);
 term(0, 4, 0, -0.08, 0.16, 0.15, 0.08, -0.04, 0.01);
 term(0, 3, 0, 0.1, -0.06, -0.07, -0.12, 0.07, -0.01);
 term(1, 4, 0, -0.31, 0.48, -0.02, 0.13, -0.03, -0.02);
 term(1, 3, 0, 0.42, -0.26, -0.38, -0.5, 0.2, -0.03);
 term(2, 4, 0, -0.7, 0.01, -0.02, -0.63, 0.0, 0.03);
 term(2, 3, 0, 2.61, -1.97, 1.74, 2.32, 0.01, 0.01);
 term(2, 2, 0, 0.32, -0.15, 0.13, 0.28, 0.0, 0.0);
 term(3, 4, 0, -0.18, 0.01, 0.0, -0.13, -0.03, 0.03);
 term(3, 3, 0, 0.75, -0.56, 0.45, 0.6, 0.08, -0.17);
 term(4, 3, 0, 0.2, -0.15, 0.1, 0.14, 0.04, -0.08);

 // perturbations by Saturn
 c[3] = cos(2*m6); s[3] = -sin(2*m6);
 term(2, 3, 0, -0.19, 0.33, 0.0, 0.0, 0.0, 0.0);

 dl = dl + (2.8 + 3.2*t);
 l = p2 * frac(0.2151379 + m1/p2 + ((5601.7+1.1*t)*t+dl)/1296.0e3);
 b = (-2522.15+(-30.18+0.04*t)*t+db) / arc;
 r = 0.3952829 +0.0000016*t + dr*1.0e-6;
 dl = 714.0+292.66*cos(m1)+71.96*cos(2*m1)+18.16*cos(3*m1)
		  +4.61*cos(4*m1)+3.81*sin(2*m1)+2.43*sin(3*m1)+1.08*sin(4*m1);
 dr = 55.94*sin(m1)+11.36*sin(2*m1)+2.6*sin(3*m1);
 db = 73.4*cos(m1)+29.82*cos(2*m1)+10.22*cos(3*m1)+3.28*cos(4*m1)
		  -40.44*sin(m1)-16.55*sin(2*m1)-5.56*sin(3*m1)-1.72*sin(4*m1);

 posvel();

return rp;
}

/*--------------------- Venus ------------------------------------------*/

Vec3 Plan200::Venus (double t)   // position of Venus at time t
 {
 /*
	 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
	 of Venus for Equinox of Date.
	 t is the time in Julian centuries since J2000.
  */

	 const double arc = 206264.81;    // 3600*180/pi = ''/r
	 const double p2 = M_PI * 2.0;

	 int i;

	 tt = t;
	 dl = 0.0; dr = 0.0; db = 0.0;
	 m1 = p2 * frac(0.4861431+415.2018375*t);
	 m2 = p2 * frac(0.1400197+162.5494552*t);
	 m3 = p2 * frac(0.9944153+99.9982208*t);
	 m4 = p2 * frac(0.0556297+53.1674631*t);
	 m5 = p2 * frac(0.0567028+8.4305083*t);
	 m6 = p2 * frac(0.8830539+3.3947206*t);

	 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m2); s3[1] = sin(m2);
	 for (i=2; i<9; i++)
		addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

	 // perturbations by Mercury
	 c[8] = 1.0; s[8] = 0.0; c[7] = cos(m1); s[7] = -sin(m1);
	 addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
	 term(1, 7, 0, 0.0, 0.0, 0.06, -0.09, 0.01, 0.0);
	 term(2, 7, 0, 0.25, -0.09, -0.09, -0.27, 0.0, 0.0);
	 term(4, 6, 0, -0.07, -0.08, -0.14, 0.14, -0.01, -0.01);
	 term(5, 6, 0, -0.35, 0.08, 0.02, 0.09, 0.0, 0.0);

	 // perturbations by Earth
	 c[7] = cos(m3); s[7] = -sin(m3);
	 for (i=7; i>0; i--)
		addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
	 term(1, 8, 0, 2.37, 2793.23, -4899.07, 0.11, 9995.27, 7027.22);
	 term(1, 8, 1, 0.1, -19.65, 34.4, 0.22, 64.95, -86.1);
	 term(1, 8, 2, 0.06, 0.04, -0.07, 0.11, -0.55, -0.07);
	 term(2, 8, 0, -170.42, 73.13, -16.59, 0.0, 67.71, 47.56);
	 term(2, 8, 1, 0.93, 2.91, 0.23, 0.0, -0.03, -0.92);
	 term(3, 8, 0, -2.31, 0.9, -0.08, 0.0, 0.04, 2.09);
	 term(1, 7, 0, -2.38, -4.27, 3.27, -1.82, 0.0, 0.0);
	 term(1, 6, 0, 0.09, 0.0, -0.08, 0.05, -0.02, -0.25);
	 term(2, 6, 0, -9.57, -5.93, 8.57, -13.83, -0.01, -0.01);
	 term(2, 5, 0, -2.47, -2.4, 0.83, -0.95, 0.16, 0.24);
	 term(3, 6, 0, -0.09, -0.05, 0.08, -0.13, -0.28, 0.12);
	 term(3, 5, 0, 7.12, 0.32, -0.62, 13.76, -0.07, 0.01);
	 term(3, 4, 0, -0.65, -0.17, 0.18, -0.73, 0.1, 0.05);
	 term(3, 3, 0, -1.08, -0.95, -0.17, 0.22, -0.03, -0.03);
	 term(4, 5, 0, 0.06, 0.0, -0.01, 0.08, 0.14, -0.18);
	 term(4, 4, 0, 0.93, -0.46, 1.06, 2.13, -0.01, 0.01);
	 term(4, 3, 0, -1.53, 0.38, -0.64, -2.54, 0.27, 0.0);
	 term(4, 2, 0, -0.17, -0.05, 0.03, -0.11, 0.02, 0.0);
	 term(5, 3, 0, 0.18, -0.28, 0.71, 0.47, -0.02, 0.04);
	 term(5, 2, 0, 0.15, -0.14, 0.3, 0.31, -0.04, 0.03);
	 term(5, 1, 0, -0.08, 0.02, -0.03, -0.11, 0.01, 0.0);
	 term(5, 0, 0, -0.23, 0.0, 0.01, -0.04, 0.0, 0.0);
	 term(6, 2, 0, 0.01, -0.14, 0.39, 0.04, 0.0, -0.01);
	 term(6, 1, 0, 0.02, -0.05, 0.12, 0.04, -0.01, 0.01);
	 term(6, 0, 0, 0.1, -0.1, 0.19, 0.19, -0.02, 0.02);
	 term(7, 1, 0, -0.03, -0.06, 0.18, -0.08, 0.0, 0.0);
	 term(8, 0, 0, -0.03, -0.02, 0.06, -0.08, 0.0, 0.0);

	 // perturbations by Mars
	 c[7] = cos(m4); s[7] = -sin(m4);
	 for (i=7; i>5; i--)
		addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
	 term(1, 5, 0, -0.65, 1.02, -0.04, -0.02, -0.02, 0.0);
	 term(2, 6, 0, -0.05, 0.04, -0.09, -0.1, 0.0, 0.0);
	 term(2, 5, 0, -0.5, 0.45, -0.79, -0.89, 0.01, 0.03);

	 // perturbations by Jupiter
	 c[7] = cos(m5); s[7] = -sin(m5);
	 for (i=7; i>5; i--)
		addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
	 term(0, 7, 0, -0.05, 1.56, 0.16, 0.04, -0.08, -0.04);
	 term(1, 7, 0, -2.62, 1.4, -2.35, -4.4, 0.02, 0.03);
	 term(1, 6, 0, -0.47, -0.08, 0.12, -0.76, 0.04, -0.18);
	 term(2, 6, 0, -0.73, -0.51, 1.27, -1.82, -0.01, 0.01);
	 term(2, 5, 0, -0.14, -0.1, 0.25, -0.34, 0.0, 0.0);
	 term(3, 5, 0, -0.01, 0.04, -0.11, -0.02, 0.0, 0.0);

	 // perturbations by Saturn
	 c[7] = cos(m6); s[7] = -sin(m6);
	 term(0, 7, 0, 0.0, 0.21, 0.0, 0.0, 0.0, -0.01);
	 term(1, 7, 0, -0.11, -0.14, 0.24, -0.2, 0.01, 0.0);

	 dl = dl+2.74*sin(p2*(0.0764+0.4174*t))+0.27*sin(p2*(0.9201+0.3307*t));
	 dl = dl+(1.9+1.8*t);

	 l = p2 * frac(0.3654783 + m2/p2 + ((5071.2+1.1*t)*t+dl)/1296.0e3);
	 b = (-67.7+(0.04+0.01*t)*t+db)/arc;
	 r = 0.7233482-0.0000002*t+dr*1.0e-6;

	 dl = 280.0+3.79*cos(m2);
	 dr = 1.37*sin(m2);
	 db = 9.54*cos(m2)-13.57*sin(m2);

	 posvel();

	return rp;
 }

/*--------------------- Mars ------------------------------------------*/

Vec3 Plan200::Mars (double t)   // position of Mars at time t
 {
 /*
	 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
	 of Mars for Equinox of Date.
	 t is the time in Julian centuries since J2000.
  */

	 const double arc = 206264.81;    // 3600*180/pi = ''/r
	 const double p2 = M_PI * 2.0;

	 int i;

	 tt = t;
	 dl = 0.0; dr = 0.0; db = 0.0;
	 m2 = p2 * frac(0.1382208+162.5482542*t);
	 m3 = p2 * frac(0.9926208+99.9970236*t);
	 m4 = p2 * frac(0.0538553+53.1662736*t);
	 m5 = p2 * frac(0.0548944+8.4290611*t);
	 m6 = p2 * frac(0.8811167+3.3935250*t);
	 c3[2] = 1.0; s3[2] = 0.0; c3[3] = cos(m4); s3[3] = sin(m4);
	 for (i=4; i<19; i++)
		addthe(c3[i-1],s3[i-1],c3[3],s3[3],c3[i],s3[i]);
	 c3[0] = c3[4];
	 s3[0] = -s3[4];
	 c3[1] = c3[3];
	 s3[1] = -s3[3];

	 // perturbations by Venus
	 c[9] = 1.0; s[9] = 0.0; c[8] = cos(m2); s[8] = -sin(m2);
	 addthe (c[8],s[8],c[8],s[8],c[7],s[7]);
	 term(2, 8, 0, -0.01, -0.03, 0.1, -0.04, 0.0, 0.0);
	 term(3, 8, 0, 0.05, 0.1, -2.08, 0.75, 0.0, 0.0);
	 term(4, 8, 0, -0.25, -0.57, -2.58, 1.18, 0.05, -0.04);
	 term(4, 7, 0, 0.02, 0.02, 0.13, -0.14, 0.0, 0.0);
	 term(5, 8, 0, 3.41, 5.38, 1.87, -1.15, 0.01, -0.01);
	 term(5, 7, 0, 0.02, 0.02, 0.11, -0.13, 0.0, 0.0);
	 term(6, 8, 0, 0.32, 0.49, -1.88, 1.21, -0.07, 0.07);
	 term(6, 7, 0, 0.03, 0.03, 0.12, -0.14, 0.0, 0.0);
	 term(7, 8, 0, 0.04, 0.06, -0.17, 0.11, -0.01, 0.01);
	 term(7, 7, 0, 0.11, 0.09, 0.35, -0.43, -0.01, 0.01);
	 term(8, 7, 0, -0.36, -0.28, -0.2, 0.25, 0.0, 0.0);
	 term(9, 7, 0, -0.03, -0.03, 0.11, -0.13, 0.0, 0.0);

	 // perturbations by Earth
	 c[8] = cos(m3); s[8] = -sin(m3);
	 for (i=8; i>0; i--)
		addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
	 term(3, 9, 0, -5.32, 38481.97, -141856.04, 0.4, -6321.67, 1876.89);
	 term(3, 9, 1, -1.12, 37.98, -138.67, -2.93, 37.28, 117.48);
	 term(3, 9, 2, -0.32, -0.03, 0.12, -1.19, 1.04, -0.4);
	 term(4, 9, 0, 28.28, 2285.8, -6608.37, 0.0, -589.35, 174.81);
	 term(4, 9, 1, 1.64, 3.37, -12.93, 0.0, 2.89, 11.1);
	 term(4, 9, 2, 0.0, 0.0, 0.0, 0.0, 0.1, -0.03);
	 term(5, 9, 0, 5.31, 189.29, -461.81, 0.0, -61.98, 18.53);
	 term(5, 9, 1, 0.31, 0.35, -1.36, 0.0, 0.25, 1.19);
	 term(6, 9, 0, 0.81, 17.96, -38.26, 0.0, -6.88, 2.08);
	 term(6, 9, 1, 0.05, 0.04, -0.15, 0.0, 0.02, 0.14);
	 term(7, 9, 0, 0.11, 1.83, -3.48, 0.0, -0.79, 0.24);
	 term(8, 9, 0, 0.02, 0.2, -0.34, 0.0, -0.09, 0.3);
	 term(1, 8, 0, 0.09, 0.06, 0.14, -0.22, 0.02, -0.02);
	 term(2, 8, 0, 0.72, 0.49, 1.55, -2.31, 0.12, -0.1);
	 term(3, 8, 0, 7.0, 4.92, 13.93, -20.48, 0.08, -0.13);
	 term(4, 8, 0, 13.08, 4.89, -4.53, 10.01, -0.05, 0.13);
	 term(4, 7, 0, 0.14, 0.05, -0.48, -2.66, 0.01, 0.14);
	 term(5, 8, 0, 1.38, 0.56, -2.0, 4.85, -0.01, 0.19);
	 term(5, 7, 0, -6.85, 2.68, 8.38, 21.42, 0.0, 0.03);
	 term(5, 6, 0, -0.08, 0.2, 1.2, 0.46, 0.0, 0.0);
	 term(6, 8, 0, 0.16, 0.07, -0.19, 0.47, -0.01, 0.05);
	 term(6, 7, 0, -4.41, 2.14, -3.33, -7.21, -0.07, -0.09);
	 term(6, 6, 0, -0.12, 0.33, 2.22, 0.72, -0.03, -0.02);
	 term(6, 5, 0, -0.04, -0.06, -0.36, 0.23, 0.0, 0.0);
	 term(7, 7, 0, -0.44, 0.21, -0.7, -1.46, -0.06, -0.07);
	 term(7, 6, 0, 0.48, -2.6, -7.25, -1.37, 0.0, 0.0);
	 term(7, 5, 0, -0.09, -0.12, -0.66, 0.5, 0.0, 0.0);
	 term(7, 4, 0, 0.03, 0.0, 0.01, -0.17, 0.0, 0.0);
	 term(8, 7, 0, -0.05, 0.03, -0.07, -0.15, -0.01, -0.01);
	 term(8, 6, 0, 0.1, -0.96, 2.36, 0.3, 0.04, 0.0);
	 term(8, 5, 0, -0.17, -0.2, -1.09, 0.94, 0.02, -0.02);
	 term(8, 4, 0, 0.05, 0.0, 0.0, -0.3, 0.0, 0.0);
	 term(9, 6, 0, 0.01, -0.1, 0.32, 0.04, 0.02, 0.0);
	 term(9, 5, 0, 0.86, 0.77, 1.86, -2.01, 0.01, -0.01);
	 term(9, 4, 0, 0.09, -0.01, -0.05, -0.44, 0.0, 0.0);
	 term(9, 3, 0, -0.01, 0.02, 0.1, 0.08, 0.0, 0.0);
	 term(10, 5, 0, 0.2, 0.16, -0.53, 0.64, -0.01, 0.02);
	 term(10, 4, 0, 0.17, -0.03, -0.14, -0.84, 0.0, 0.01);
	 term(10, 3, 0, -0.02, 0.03, 0.16, 0.09, 0.0, 0.0);
	 term(11, 4, 0, -0.55, 0.15, 0.3, 1.1, 0.0, 0.0);
	 term(11, 3, 0, -0.02, 0.04, 0.2, 0.1, 0.0, 0.0);
	 term(12, 4, 0, -0.09, 0.03, -0.1, -0.33, 0.0, -0.01);
	 term(12, 3, 0, -0.05, 0.11, 0.48, 0.21, -0.01, 0.0);
	 term(13, 3, 0, 0.1, -0.35, -0.52, -0.15, 0.0, 0.0);
	 term(13, 2, 0, -0.01, -0.02, -0.1, 0.07, 0.0, 0.0);
	 term(14, 3, 0, 0.01, -0.04, 0.18, 0.4, 0.01, 0.0);
	 term(14, 2, 0, -0.05, -0.07, -0.29, 0.2, 0.01, 0.0);
	 term(15, 2, 0, 0.23, 0.27, 0.25, -0.21, 0.0, 0.0);
	 term(16, 2, 0, 0.02, 0.03, -0.1, 0.09, 0.0, 0.0);
	 term(16, 1, 0, 0.05, 0.01, 0.03, -0.23, 0.0, 0.3);
	 term(17, 1, 0, -1.53, 0.27, 0.06, 0.42, 0.0, 0.0);
	 term(18, 1, 0, -0.14, 0.02, -0.1, -0.55, -0.01, -0.02);
	 term(18, 0, 0, 0.03, -0.06, -0.25, -0.11, 0.0, 0.0);

	 // perturbation by Jupiter
	 c[8] = cos(m5); s[8] = -sin(m5);
	 for (i=8; i>4; i--)
		addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
	 term(0, 8, 0, 0.05, 0.03, 0.08, -0.14, 0.01, -0.01);
	 term(1, 8, 0, 0.39, 0.27, 0.92, -1.5, -0.03, -0.06);
	 term(1, 7, 0, -0.16, 0.03, 0.13, 0.67, -0.01, 0.06);
	 term(1, 6, 0, -0.02, 0.01, 0.05, 0.09, 0.0, 0.01);
	 term(2, 8, 0, 3.56, 1.13, -5.41, -7.18, -0.25, -0.24);
	 term(2, 7, 0, -1.44, 0.25, 1.24, 7.96, 0.02, 0.31);
	 term(2, 6, 0, -0.21, 0.11, 0.55, 1.04, 0.01, 0.05);
	 term(2, 5, 0, -0.02, 0.02, 0.11, 0.11, 0.0, 0.01);
	 term(3, 8, 0, 16.67, -19.15, 61.0, 53.36, -0.06, -0.07);
	 term(3, 7, 0, -21.64, 3.18, -7.77, -54.64, -0.31, 0.5);
	 term(3, 6, 0, -2.82, 1.45, -2.53, -5.73, 0.01, 0.07);
	 term(3, 5, 0, -0.31, 0.28, -0.34, -0.51, 0.0, 0.0);
	 term(4, 8, 0, 2.15, -2.29, 7.04, 6.94, 0.33, 0.19);
	 term(4, 7, 0, -15.69, 3.31, -15.7, -73.17, -0.17, -0.25);
	 term(4, 6, 0, -1.73, 1.95, -9.19, -7.2, 0.02, -0.03);
	 term(4, 5, 0, -0.01, 0.33, -1.42, 0.08, 0.01, -0.01);
	 term(4, 4, 0, 0.03, 0.03, -0.13, 0.12, 0.0, 0.0);
	 term(5, 8, 0, 0.26, -0.28, 0.73, 0.71, 0.08, 0.04);
	 term(5, 7, 0, -2.06, 0.46, -1.61, -6.72, -0.13, -0.25);
	 term(5, 6, 0, -1.28, -0.27, 2.21, -6.9, -0.04, -0.02);
	 term(5, 5, 0, -0.22, 0.08, -0.44, -1.25, 0.0, 0.01);
	 term(5, 4, 0, -0.02, 0.03, -0.15, -0.08, 0.0, 0.0);
	 term(6, 8, 0, 0.03, -0.03, 0.08, 0.08, 0.01, 0.01);
	 term(6, 7, 0, -0.26, 0.06, -0.17, -0.7, -0.03, -0.05);
	 term(6, 6, 0, -0.2, -0.05, 0.22, -0.79, -0.01, -0.02);
	 term(6, 5, 0, -0.11, -0.14, 0.93, -0.6, 0.0, 0.0);
	 term(6, 4, 0, -0.04, -0.02, 0.09, -0.23, 0.0, 0.0);
	 term(7, 5, 0, -0.02, -0.03, 0.13, -0.09, 0.0, 0.0);
	 term(7, 4, 0, 0.0, -0.03, 0.21, 0.01, 0.0, 0.0);

	 // perturbations by Saturn
	 c[8] = cos(m6); s[8] = -sin(m6);
	 for (i=8; i>5; i--)
		addthe(c[i],s[i],c[8],s[8],c[i-1],s[i-1]);
	 term(1, 8, 0, 0.03, 0.13, 0.48, -0.13, 0.02, 0.0);
	 term(2, 8, 0, 0.27, 0.84, 0.4, -0.43, 0.01, -0.01);
	 term(2, 7, 0, 0.12, -0.04, -0.33, -0.55, -0.01, -0.02);
	 term(2, 6, 0, 0.02, -0.01, -0.07, -0.08, 0.0, 0.0);
	 term(3, 8, 0, 1.12, 0.76, -2.66, 3.91, -0.01, 0.01);
	 term(3, 7, 0, 1.49, -0.95, 3.07, 4.83, 0.04, -0.05);
	 term(3, 6, 0, 0.21, -0.18, 0.55, 0.64, 0.0, 0.0);
	 term(4, 8, 0, 0.12, 0.1, -0.29, 0.34, -0.01, 0.02);
	 term(4, 7, 0, 0.51, -0.36, 1.61, 2.25, 0.03, 0.01);
	 term(4, 6, 0, 0.1, -0.1, 0.5, 0.43, 0.0, 0.0);
	 term(4, 5, 0, 0.01, -0.02, 0.11,  0.05, 0.0, 0.0);
	 term(5, 7, 0, 0.07, -0.05, 0.16, 0.22, 0.01, 0.01);

	 dl = dl+52.49*sin(p2*(0.1868+0.0549*t))+0.61*sin(p2*(0.922+0.3307*t))
			  +0.32*sin(p2*(0.4731+2.1485*t))+0.28*sin(p2*(0.9467+1.1133*t));
	 dl = dl + (0.14+0.87*t-0.11*t*t);
	 l = p2 * frac(0.9334591+m4/p2+((6615.5+1.1*t)*t+dl)/1296.0e3);
	 r = 1.5303352 + 0.0000131*t + dr*1.0e-6;
	 b = (596.32+(-2.92-0.1*t)*t+db) / arc;

	 dl = 91.50 + 17.07*cos(m4) + 2.03*cos(2*m4);
	 dr = 12.98*sin(m4) + 1.21*cos(2*m4);
	 db = 0.83*cos(m4) + 2.8*sin(m4);

	 posvel();

	return rp;
 }

/*--------------------- Jupiter ------------------------------------------*/

Vec3 Plan200::Jupiter (double t)   // position of Jupiter at time t
 {
 /*
	 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
	 of Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
  */

	 const double arc = 206264.81;    // 3600*180/pi = ''/r
	 const double p2 = M_PI * 2.0;

	 int i;

	 tt = t;
	 dl = 0.0; dr = 0.0; db = 0.0;
	 m5 = p2 * frac(0.0565314+8.4302963*t);
	 m6 = p2 * frac(0.8829867+3.3947688*t);
	 m7 = p2 * frac(0.3969537+1.1902586*t);
	 c3[1] = 1.0; s3[1] = 0.0;
	 c3[2] = cos(m5); s3[2] = sin(m5); c3[0] = c3[2]; s3[0] = -s3[2];
	 for (i=3; i<7; i++)
		addthe(c3[i-1],s3[i-1],c3[2],s3[2],c3[i],s3[i]);

	 // perturbations by Saturn
	 c[10] = 1.0; s[10] = 0.0; c[9] = cos(m6); s[9] = -sin(m6);
	 for (i=9; i>0; i--)
		addthe(c[i],s[i],c[9],s[9],c[i-1],s[i-1]);
	 term(0, 9, 0, -0.2, 1.4, 2.0, 0.6, 0.1, -0.2);
	 term(1, 9, 0, 9.4, 8.9, 3.9, -8.3, -0.4, -1.4);
	 term(1, 8, 0, 5.6, -3.0, -5.4, -5.7, -2.0, 0.0);
	 term(1, 7, 0, -4.0, -0.1, 0.0, 5.5, 0.0, 0.0);
	 term(1, 5, 0, 3.3, -1.6, -1.6, -3.1, -0.5, -1.2);
	 term(2, 10, 0, -113.1, 19998.6, -25208.2, -142.2, -4670.7, 288.9);
	 term(2, 10, 1, -76.1, 66.9, -84.2, -95.8, 21.6, 29.4);
	 term(2, 10, 2, -0.5, -0.3, 0.4, -0.7, 0.1, -0.1);
	 term(2, 9, 0, 78.8, -14.5, 11.5, 64.4, -0.2, 0.2);
	 term(2, 8, 0, -2.0, -132.4, 28.8, 4.3, -1.7, 0.4);
	 term(2, 8, 1, -1.1, -0.7, 0.2, -0.3, 0.0, 0.0);
	 term(2, 7, 0, -7.5, -6.8, -0.4, -1.1, 0.6, -0.9);
	 term(2, 6, 0, 0.7, 0.7, 0.6, -1.1, 0.0, -0.2);
	 term(2, 5, 0, 51.5, -26.0, -32.5, -64.4, -4.9, -12.4);
	 term(2, 5, 1, -1.2, -2.2, -2.7, 1.5, -0.4, 0.3);
	 term(3, 10, 0, -3.4, 632.0, -610.6, -6.5, -226.8, 12.7);
	 term(3, 10, 1, -4.2, 3.8, -4.1, -4.5, 0.2, 0.6);
	 term(3, 9, 0, 5.3, -0.7, 0.7, 6.1, 0.2, 1.1);
	 term(3, 8, 0, -76.4, -185.1, 260.2, -108.0, 1.6, 0.0);
	 term(3, 7, 0, 66.7, 47.8, -51.4, 69.8, 0.9, 0.3);
	 term(3, 7, 1, 0.6, -1.0, 1.0, 0.6, 0.0, 0.0);
	 term(3, 6, 0, 17.0, 1.4, -1.8, 9.6, 0.0, -0.1);
	 term(3, 5, 0, 1066.2, -518.3, -1.3, -23.9, 1.8, -0.3);
	 term(3, 5, 1, -25.4, -40.3, -0.9, 0.3, 0.0, 0.0);
	 term(3, 5, 2, -0.7, 0.5, 0.0, 0.0, 0.0, 0.0);
	 term(4, 10, 0, -0.1, 28.0, -22.1, -0.2, -12.5, 0.7);
	 term(4, 8, 0, -5.0, -11.5, 11.7, -5.4, 2.1, -1.0);
	 term(4, 7, 0, 16.9, -6.4, 13.4, 26.9, -0.5, 0.8);
	 term(4, 6, 0, 7.2, -13.3, 20.9, 10.5, 0.1, -0.1);
	 term(4, 5, 0, 68.5, 134.3, -166.9, 86.5, 7.1, 15.2);
	 term(4, 5, 1, 3.5, -2.7, 3.4, 4.3, 0.5, -0.4);
	 term(4, 4, 0, 0.6, 1.0, -0.9, 0.5, 0.0, 0.0);
	 term(4, 3, 0, -1.1, 1.7, -0.4, -0.2, 0.0, 0.0);
	 term(5, 10, 0, 0.0, 1.4, -1.0, 0.0, -0.6, 0.0);
	 term(5, 8, 0, -0.3, -0.7, 0.4, -0.2, 0.2, -0.1);
	 term(5, 7, 0, 1.1, -0.6, 0.9, 1.2, 0.1, 0.2);
	 term(5, 6, 0, 3.2, 1.7, -4.1, 5.8, 0.2, 0.1);
	 term(5, 5, 0, 6.7, 8.7, -9.3, 8.7, -1.1, 1.6);
	 term(5, 4, 0, 1.5, -0.3, 0.6, 2.4, 0.0, 0.0);
	 term(5, 3, 0, -1.9, 2.3, -3.2, -2.7, 0.0, -0.1);
	 term(5, 2, 0, 0.4, -1.8, 1.9, 0.5, 0.0, 0.0);
	 term(5, 1, 0, -0.2, -0.5, 0.3, -0.1, 0.0, 0.0);
	 term(5, 0, 0, -8.6, -6.8, -0.4, 0.1, 0.0, 0.0);
	 term(5, 0, 1, -0.5, 0.6, 0.0, 0.0, 0.0, 0.0);
	 term(6, 5, 0, -0.1, 1.5, -2.5, -0.8, -0.1, 0.1);
	 term(6, 4, 0, 0.1, 0.8, -1.6, 0.1, 0.0, 0.0);
	 term(6, 1, 0, -0.5, -0.1, 0.1, -0.8, 0.0, 0.0);
	 term(6, 0, 0, 2.5, -2.2, 2.8, 3.1, 0.1, -0.2);

	 // perturbations by Uranus
	 c[9] = cos(m7); s[9] = -sin(m7);
	 addthe(c[9],s[9],c[9],s[9],c[8],s[8]);
	 term(2, 9, 0, 0.4, 0.9, 0.0, 0.0, 0.0, 0.0);
	 term(2, 8, 0, 0.4, 0.4, -0.4, 0.3, 0.0, 0.0);

	 // perturbations by Saturn and Uranus
	 double phi, x, y;

	 phi = (2*m5-6*m6+3*m7); x = cos(phi); y = sin(phi);
	 dl = dl-0.8*x+8.5*y; dr = dr-0.1*x;
	 addthe(x,y,c3[2],s3[2],x,y);
	 dl = dl+0.4*x+0.5*y; dr = dr-0.7*x+0.5*y; db = db-0.1*x;

	 l = p2 * frac(0.0388910 + m5/p2 + ((5025.2+0.8*t)*t+dl)/1296.0e3);
	 b = (227.3 - 0.3*t + db) /  arc;
	 r = 5.208873 + 0.000041*t + dr*1.0e-5;

	 dl = 14.5+1.41*cos(m5); dr = 3.66*sin(m5); db = 0.33*sin(m5);

	 posvel();

	return rp;
 }


/*--------------------- Saturn ------------------------------------------*/

Vec3 Plan200::Saturn (double t)   // position of Saturn at time t
 {
 /*
	 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
	 of Saturn for Equinox of Date.
	 t is the time in Julian centuries since J2000.
  */

	 const double arc = 206264.81;    // 3600*180/pi = ''/r
	 const double p2 = M_PI * 2.0;

	 int i;

	 tt = t;
	 dl = 0.0; dr = 0.0; db = 0.0;
	 m5 = p2 * frac(0.0565314+8.4302963*t);
	 m6 = p2 * frac(0.8829867+3.3947688*t);
	 m7 = p2 * frac(0.3969537+1.1902586*t);
	 m8 = p2 * frac(0.7208473+0.6068623*t);
	 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m6); s3[1] = sin(m6);
	 for (i=2; i<12; i++)
		addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

	 // perturbations by Jupiter
	 c[6] = 1.0; s[6] = 0.0; c[7] = cos(m5); s[7] = sin(m5);
	 for (i=6; i>0; i--)
		addthe(c[i],s[i],c[7],-s[7],c[i-1],s[i-1]);
	 term(0, 5, 0, 12.0, -1.4, -13.9, 6.4, 1.2, -1.8);
	 term(0, 4, 0, 0.0, -0.2, -0.9, 1.0, 0.0, -0.1);
	 term(1, 7, 0, 0.9, 0.4, -1.8, 1.9, 0.2, 0.2);
	 term(1, 6, 0, -348.3, 22907.7, -52915.5, -752.2, -3266.5, 8314.4);
	 term(1, 6, 1, -225.2, -146.2, 337.7, -521.3, 79.6, 17.4);
	 term(1, 6, 2, 1.3, -1.4, 3.2, 2.9, 0.1, -0.4);
	 term(1, 5, 0, -1.0, -30.7, 108.6, -815.0, -3.6, -9.3);
	 term(1, 4, 0, -2.0, -2.7, -2.1, -11.9, -0.1, -0.4);
	 term(2, 7, 0, 0.1, 0.2, -1.0, 0.3, 0.0, 0.0);
	 term(2, 6, 0, 44.2, 724.0, -1464.3, -34.7, -188.7, 459.1);
	 term(2, 6, 1, -17.0, -11.3, 18.9, -28.9, 1.0, -3.7);
	 term(2, 5, 0, -3.5, -426.6, -546.5, -26.5, -1.6, -2.7);
	 term(2, 5, 1, 3.5, -2.2, -2.6, -4.3, 0.0, 0.0);
	 term(2, 4, 0, 10.5, -30.9, -130.5, -52.3, -1.9, 0.2);
	 term(2, 3, 0, -0.2, -0.4, -1.2, -0.1, -0.1, 0.0);
	 term(3, 6, 0, 6.5, 30.5, -61.1, 0.4, -11.6, 28.1);
	 term(3, 6, 1, -1.2, -0.7, 1.1, -1.8, -0.2, -0.6);
	 term(3, 5, 0, 29.0, -40.2, 98.2, 45.3, 3.2, -9.4);
	 term(3, 5, 1, 0.6, 0.6, -1.0, 1.3, 0.0, 0.0);
	 term(3, 4, 0, -27.0, -21.1, -68.5, 8.1, -19.8, 5.4);
	 term(3, 4, 1, 0.9, -0.5, -0.4, -2.0, -0.1, -0.8);
	 term(3, 3, 0, -5.4, -4.1, -19.1, 26.2, -0.1, -0.1);
	 term(4, 6, 0, 0.6, 1.4, -3.0, -0.2, -0.6, 1.6);
	 term(4, 5, 0, 1.5, -2.5, 12.4, 4.7, 1.0, -1.1);
	 term(4, 4, 0, -821.9, -9.6, -26.0, 1873.6, -70.5, -4.4);
	 term(4, 4, 1, 4.1, -21.9, -50.3, -9.9, 0.7, -3.0);
	 term(4, 3, 0, -2.0, -4.7, -19.3, 8.2, -0.1, -0.3);
	 term(4, 2, 0, -1.5, 1.3, 6.5, 7.3, 0.0, 0.0);
	 term(5, 4, 0, -2627.6, -1277.3, 117.4, -344.1, -13.8, -4.3);
	 term(5, 4, 1, 63.0, -98.6, 12.7, 6.7, 0.1, -0.2);
	 term(5, 4, 2, 1.7, 1.2, -0.2, 0.3, 0.0, 0.0);
	 term(5, 3, 0, 0.4, -3.6, -11.3, -1.6, 0.0, -0.3);
	 term(5, 2, 0, -1.4, 0.3, 1.5, 6.3, -0.1, 0.0);
	 term(5, 1, 0, 0.3, 0.6, 3.0, -1.7, 0.0, 0.0);
	 term(6, 4, 0, -146.7, -73.7, 166.4, -334.3, -43.6, -46.7);
	 term(6, 4, 1, 5.2, -6.8, 15.1, 11.4, 1.7, -1.0);
	 term(6, 3, 0, 1.5, -2.9, -2.2, -1.3, 0.1, -0.1);
	 term(6, 2, 0, -0.7, -0.2, -0.7, 2.8, 0.0, 0.0);
	 term(6, 1, 0, 0.0, 0.5, 2.5, -0.1, 0.0, 0.0);
	 term(6, 0, 0, 0.3, -0.1, -0.3, -1.2, 0.0, 0.0);
	 term(7, 4, 0, -9.6, -3.9, 9.6, -18.6, -4.7, -5.3);
	 term(7, 4, 1, 0.4, -0.5, 1.0, 0.9, 0.3, -0.1);
	 term(7, 3, 0, 3.0, 5.3, 7.5, -3.5, 0.0, 0.0);
	 term(7, 2, 0, 0.2, 0.4, 1.6, -1.3, 0.0, 0.0);
	 term(7, 1, 0, -0.1, 0.2, 1.0, 0.5, 0.0, 0.0);
	 term(7, 0, 0, 0.2, 0.0, 0.2, -1.0, 0.0, 0.0);
	 term(8, 4, 0, -0.7, -0.2, 0.6, -1.2, -0.4, -0.4);
	 term(8, 3, 0, 0.5, 1.0, -2.0, 1.5, 0.1, 0.2);
	 term(8, 2, 0, 0.4, 1.3, 3.6, -0.9, 0.0, -0.1);
	 term(9, 2, 0, 4.0, -8.7, -19.9, -9.9, 0.2, -0.4);
	 term(9, 2, 1, 0.5, 0.3, 0.8, -1.8, 0.0, 0.0);
	 term(10, 2, 0, 21.3, -16.8, 3.3, 3.3, 0.2, -0.2);
	 term(10, 2, 1, 1.0, 1.7, -0.4, 0.4, 0.0, 0.0);
	 term(11, 2, 0, 1.6, -1.3, 3.0, 3.7, 0.8, -0.2);

	 // perturbations by Uranus
	 c[5] = cos(m7); s[5] = -sin(m7);
	 for (i=5; i>1; i--)
		addthe(c[i],s[i],c[5],s[5],c[i-1],s[i-1]);
	 term(0, 5, 0, 1.0, 0.7, 0.4, -1.5, 0.1, 0.0);
	 term(0, 4, 0, 0.0, -0.4, -1.1, 0.1, -0.1, -0.1);
	 term(0, 3, 0, -0.9, -1.2, -2.7, 2.1, -0.5, -0.3);
	 term(1, 5, 0, 7.8, -1.5, 2.3, 12.7, 0.0, 0.0);
	 term(1, 4, 0, -1.1, -8.1, 5.2, -0.3, -0.3, -0.3);
	 term(1, 3, 0, -16.4, -21.0, -2.1, 0.0, 0.4, 0.0);
	 term(2, 5, 0, 0.6, -0.1, 0.1, 1.2, 0.1, 0.0);
	 term(2, 4, 0, -4.9, -11.7, 31.5, -13.3, 0.0, -0.2);
	 term(2, 3, 0, 19.1, 10.0, -22.1, 42.1, 0.1, -1.1);
	 term(2, 2, 0, 0.9, -0.1, 0.1, 1.4, 0.0, 0.0);
	 term(3, 4, 0, -0.4, -0.9, 1.7, -0.8, 0.0, -0.3);
	 term(3, 3, 0, 2.3, 0.0, 1.0, 5.7, 0.3, 0.3);
	 term(3, 2, 0, 0.3, -0.7, 2.0, 0.7, 0.0, 0.0);
	 term(3, 1, 0, -0.1, -0.4, 1.1, -0.3, 0.0, 0.0);

	 // perturbations by Neptune
	 c[5] = cos(m8); s[5] = -sin(m8);
	 addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
	 term(1, 5, 0, -1.3, -1.2, 2.3, -2.5, 0.0, 0.0);
	 term(1, 4, 0, 1.0, -0.1, 0.1, 1.4, 0.0, 0.0);
	 term(2, 4, 0, 1.1, -0.1, 0.2, 3.3, 0.0, 0.0);

	 // perturbations by Jupiter and Uranus
	 double phi, x, y;

	 phi = (-2*m5+5*m6-3*m7); x = cos(phi); y = sin(phi);
	 dl = dl-0.8*x-0.1*y; dr = dr-0.2*x+1.8*y; db = db+0.3*x+0.5*y;
	 addthe(x,y,c3[1],s3[1],x,y);
	 dl = dl+(2.4-0.7*t)*x+(27.8-0.4*t)*y;
	 dr = dr+2.1*x-0.2*y;
	 addthe(x,y,c3[1],s3[1],x,y);
	 dl = dl+0.1*x+1.6*y; dr = dr-3.6*x+0.3*y; db = db-0.2*x+0.6*y;

	 l = p2 * frac(0.2561136+m6/p2+((5018.6+t*1.9)*t+dl)/1296.0e3);
	 b = (175.1 - 10.2*t + db) / arc;
	 r = 9.557584 - 0.000186*t + dr*1.0e-5;

	 dl = 5.84+0.65*cos(m6); dr = 3.09*sin(m6); db = 0.24*cos(m6);

	 posvel();

	return rp;
 }


/*--------------------- Uranus ------------------------------------------*/

Vec3 Plan200::Uranus (double t)   // position of Uranus at time t
 {
 /*
	 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
	 of Uranus for Equinox of Date.
	 t is the time in Julian centuries since J2000.
  */

	 const double arc = 206264.81;    // 3600*180/pi = ''/r
	 const double p2 = M_PI * 2.0;

	 int i;

	 tt = t;
	 dl = 0.0; dr = 0.0; db = 0.0;
	 m5 = p2 * frac(0.0564472+8.4302889*t);
	 m6 = p2 * frac(0.8829611+3.3947583*t);
	 m7 = p2 * frac(0.3967117+1.1902849*t);
	 m8 = p2 * frac(0.7216833+0.6068528*t);
	 c3[2] = 1.0; s3[2] = 0.0; c3[3] = cos(m7); s3[3] = sin(m7);
	 for (i=4; i<10; i++)
		addthe(c3[i-1],s3[i-1],c3[3],s3[3],c3[i],s3[i]);
	 c3[1] = c3[3];
	 s3[1] = -s3[3];
	 c3[0] = c3[4];
	 s3[0] = -s3[4];

	 // perturbations by Jupiter
	 c[8] = 1.0; s[8] = 0.0; c[7] = cos(m5); s[7] = -sin(m5);
	 addthe(c[7],s[7],c[7],s[7],c[6],s[6]);
	 term(1, 7, 0, 0.0, 0.0, -0.1, 1.7, -0.1, 0.0);
	 term(2, 7, 0, 0.5, -1.2, 18.9, 9.1, -0.9, 0.1);
	 term(3, 7, 0, -21.2, 48.7, -455.5, -198.8, 0.0, 0.0);
	 term(3, 6, 0, -0.5, 1.2, -10.9, -4.8, 0.0, 0.0);
	 term(4, 7, 0, -1.3, 3.2, -23.2, -11.1, 0.3, 0.1);
	 term(4, 6, 0, -0.2, 0.2, 1.1, 1.5, 0.0, 0.0);
	 term(5, 7, 0, 0.0, 0.2, -1.8, 0.4, 0.0, 0.0);

	 // perturbations by Saturn
	 c[7] = cos(m6); s[7] = -sin(m6);
	 for (i=7; i>4; i--)
		addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
	 term(2, 7, 0, 1.4, -0.5, -6.4, 9.0, -0.4, -0.8);
	 term(3, 7, 0, -18.6, -12.6, 36.7, -336.8, 1.0, 0.3);
	 term(3, 6, 0, -0.7, -0.3, 0.5, -7.5, 0.1, 0.0);
	 term(4, 7, 0, 20.0, -141.6, -587.1, -107.0, 3.1, -0.8);
	 term(4, 7, 1, 1.0, 1.4, 5.8, -4.0, 0.0, 0.0);
	 term(4, 6, 0, 1.6, -3.8, -35.6, -16.0, 0.0, 0.0);
	 term(5, 7, 0, 75.3, -100.9, 128.9, 77.5, -0.8, 0.1);
	 term(5, 7, 1, 0.2, 1.8, -1.9, 0.3, 0.0, 0.0);
	 term(5, 6, 0, 2.3, -1.3, -9.5, -17.9, 0.0, 0.1);
	 term(5, 5, 0, -0.7, -0.5, -4.9, 6.8, 0.0, 0.0);
	 term(6, 7, 0, 3.4, -5.0, 21.6, 14.3, -0.8, -0.5);
	 term(6, 6, 0, 1.9, 0.1, 1.2, -12.1, 0.0, 0.0);
	 term(6, 5, 0, -0.1, -0.4, -3.9, 1.2, 0.0, 0.0);
	 term(6, 4, 0, -0.2, 0.1, 1.6, 1.8, 0.0, 0.0);
	 term(7, 7, 0, 0.2, -0.3, 1.0, 0.6, -0.1, 0.0);
	 term(7, 6, 0, -2.2, -2.2, -7.7, 8.5, 0.0, 0.0);
	 term(7, 5, 0, 0.1, -0.2, -1.4, -0.4, 0.0, 0.0);
	 term(7, 4, 0, -0.1, 0.0, 0.1, 1.2, 0.0, 0.0);
	 term(8, 6, 0, -0.2, -0.6, 1.4, -0.7, 0.0, 0.0);

	 // perturbations by Neptune
	 c[7] = cos(m8); s[7] = -sin(m8);
	 for (i=7; i>0; i--)
		addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
	 term(3, 8, 0, -78.1, 19518.1, -90718.2, -334.7, 2759.5, -311.9);
	 term(3, 8, 1, -81.6, 107.7, -497.4, -379.5, -2.8, -43.7);
	 term(3, 8, 2, -6.6, -3.1, 14.4, -30.6, -0.4, -0.5);
	 term(3, 8, 3, 0.0, -0.5, 2.4, 0.0, 0.0, 0.0);
	 term(4, 8, 0, -2.4, 586.1, -2145.2, -15.3, 130.6, -14.3);
	 term(4, 8, 1, -4.5, 6.6, -24.2, -17.8, 0.7, -1.6);
	 term(4, 8, 2, -0.4, 0.0, 0.1, -1.4, 0.0, 0.0);
	 term(5, 8, 0, 0.0, 24.5, -76.2, -0.6, 7.0, -0.7);
	 term(5, 8, 1, -0.2, 0.4, -1.4, -0.8, 0.1, -0.1);
	 term(6, 8, 0, 0.0, 1.1, -3.0, 0.1, 0.4, 0.0);
	 term(1, 7, 0, -0.2, 0.2, 0.7, 0.7, -0.1, 0.0);
	 term(2, 7, 0, -2.8, 2.5, 8.7, 10.5, -0.4, -0.1);
	 term(3, 7, 0, -28.4, 20.3, -51.4, -72.0, 0.0, 0.0);
	 term(3, 6, 0, -0.6, -0.1, 4.2, -14.6, 0.2, 0.4);
	 term(3, 5, 0, 0.2, 0.5, 3.4, -1.6, -0.1, 0.1);
	 term(4, 7, 0, -1.8, 1.3, -5.5, -7.7, 0.0, 0.3);
	 term(4, 6, 0, 29.4, 10.2, -29.0, 83.2, 0.0, 0.0);
	 term(4, 5, 0, 8.8, 17.8, -41.9, 21.5, -0.1, -0.3);
	 term(4, 4, 0, 0.0, 0.1, -2.1, -0.9, 0.1, 0.0);
	 term(5, 6, 0, 1.5, 0.5, -1.7, 5.1, 0.1, -0.2);
	 term(5, 5, 0, 4.4, 14.6, -84.3, 25.2, 0.1, -0.1);
	 term(5, 4, 0, 2.4, -4.5, 12.0, 6.2, 0.0, 0.0);
	 term(5, 3, 0, 2.9, -0.9, 2.1, 6.2, 0.0, 0.0);
	 term(6, 5, 0, 0.3, 1.0, -4.0, 1.1, 0.1, -0.1);
	 term(6, 4, 0, 2.1, -2.7, 17.9, 14.0, 0.0, 0.0);
	 term(6, 3, 0, 3.0, -0.4, 2.3, 17.6, -0.1, -0.1);
	 term(6, 2, 0, -0.6, -0.5, 1.1, -1.6, 0.0, 0.0);
	 term(7, 4, 0, 0.2, -0.2, 1.0, 0.8, 0.0, 0.0);
	 term(7, 3, 0, -0.9, -0.1, 0.6, -7.1, 0.0, 0.0);
	 term(7, 2, 0, -0.5, -0.6, 3.8, -3.6, 0.0, 0.0);
	 term(7, 1, 0, 0.0, -0.5, 3.0, 0.1, 0.0, 0.0);
	 term(8, 2, 0, 0.2, 0.3, -2.7, 1.6, 0.0, 0.0);
	 term(8, 1, 0, -0.1, 0.2, -2.0, -0.4, 0.0, 0.0);
	 term(9, 1, 0, 0.1, -0.2, 1.3, 0.5, 0.0, 0.0);
	 term(9, 0, 0, 0.1, 0.0, 0.4, 0.9, 0.0, 0.0);

	 // perturbations by Jupiter and Saturn
	 c[7] = cos(m6); s[7] = -sin(m6);
	 c[4] = cos(-4*m6+2*m5); s[4] = sin(-4*m6+2*m5);
	 for (i=4; i>2; i--)
		addthe(c[i],s[i],c[7],s[7],c[i-1],s[i-1]);
	 term(0, 4, 0, -0.7, 0.4, -1.5, -2.5, 0.0, 0.0);
	 term(1, 4, 0, -0.1, -0.1, -2.2, 1.0, 0.0, 0.0);
	 term(3, 3, 0, 0.1, -0.4, 1.4, 0.2, 0.0, 0.0);
	 term(3, 2, 0, 0.4, 0.5, -0.8, -0.8, 0.0, 0.0);
	 term(4, 2, 0, 5.7, 6.3, 28.5, -25.5, 0.0, 0.0);
	 term(4, 2, 1, 0.1, -0.2, -1.1, -0.6, 0.0, 0.0);
	 term(5, 2, 0, -1.4, 29.2, -11.4, 1.1, 0.0, 0.0);
	 term(5, 2, 1, 0.8, -0.4, 0.2, 0.3, 0.0, 0.0);
	 term(6, 2, 0, 0.0, 1.3, -6.0, -0.1, 0.0, 0.0);

	 l = p2 * frac(0.4734843+m7/p2+((5082.3+34.2*t)*t+dl)/1296.0E3);
	 b = (-130.61+(-0.54+0.04*t)*t+db)/arc;
	 r = 19.211991+(-0.000333-0.000005*t)*t+dr*1.0e-5;

	 dl = 2.05+0.19*cos(m7); dr = 1.86*sin(m7); db = -0.03*sin(m7);

	 posvel();

	return rp;
 }


/*--------------------- Neptune ------------------------------------------*/

Vec3 Plan200::Neptune (double t)   // position of Neptune at time t
 {
 /*
	 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
	 of Neptune for Equinox of Date.
	 t is the time in Julian centuries since J2000.
  */

	 const double arc = 206264.81;    // 3600*180/pi = ''/r
	 const double p2 = M_PI * 2.0;

	 int i;

	 tt = t;
	 dl = 0.0; dr = 0.0; db = 0.0;
	 m5 = p2 * frac(0.0563867+8.4298907*t);
	 m6 = p2 * frac(0.8825086+3.3957748*t);
	 m7 = p2 * frac(0.3965358+1.1902851*t);
	 m8 = p2 * frac(0.7214906+0.6068526*t);
	 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m8); s3[1] = sin(m8);
	 for (i=2; i<7; i++)
		addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

	 // perturbations by Jupiter
	 c[6] = 1.0; s[6] = 0.0; c[5] = cos(m5); s[5] = -sin(m5);
	 addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
	 term(0, 5, 0, 0.1, 0.1, -3.0, 1.8, -0.3, -0.3);
	 term(1, 6, 0, 0.0, 0.0, -15.9, 9.0, 0.0, 0.0);
	 term(1, 5, 0, -17.6, -29.3, 416.1, -250.0, 0.0, 0.0);
	 term(1, 4, 0, -0.4, -0.7, 10.4, -6.2, 0.0, 0.0);
	 term(2, 5, 0, -0.2, -0.4, 2.4, -1.4, 0.4, -0.3);

	 // perturbations by Saturn
	 c[6] = 1.0; s[6] = 0.0; c[5] = cos(m6); s[5] = -sin(m6);
	 addthe(c[5],s[5],c[5],s[5],c[4],s[4]);
	 term(0, 5, 0, -0.1, 0.0, 0.2, -1.8, -0.1, -0.5);
	 term(1, 6, 0, 0.0, 0.0, -8.3, -10.4, 0.0, 0.0);
	 term(1, 5, 0, 13.6, -12.7, 187.5, 201.1, 0.0, 0.0);
	 term(1, 4, 0, 0.4, -0.4, 4.5, 4.5, 0.0, 0.0);
	 term(2, 5, 0, 0.4, -0.1, 1.7, -3.2, 0.2, 0.2);
	 term(2, 4, 0, -0.1, 0.0, -0.2, 2.7, 0.0, 0.0);

	 // perturbations by Uranus
	 c[6] = 1.0; s[6] = 0.0; c[5] = cos(m7); s[5] = -sin(m7);
	 for (i=5; i>0; i--)
		addthe(c[i],s[i],c[5],s[5],c[i-1],s[i-1]);
	 term(1, 6, 0, 32.3, 3549.5, -25880.2, 235.8, -6360.5, 374.0);
	 term(1, 6, 1, 31.2, 34.4, -251.4, 227.4, 34.9, 29.3);
	 term(1, 6, 2, -1.4, 3.9, -28.6, -10.1, 0.0, -0.9);
	 term(2, 6, 0, 6.1, 68.0, -111.4, 2.0, -54.7, 3.7);
	 term(2, 6, 1, 0.8, -0.2, -2.1, 2.0, -0.2, 0.8);
	 term(3, 6, 0, 0.1, 1.0, -0.7, 0.0, -0.8, 0.1);
	 term(0, 5, 0, -0.1, -0.3, -3.6, 0.0, 0.0, 0.0);
	 term(1, 6, 0, 0.0, 0.0, 5.5, -6.9, 0.1, 0.0);
	 term(1, 5, 0, -2.2, -1.6, -116.3, 163.6, 0.0, -0.1);
	 term(1, 4, 0, 0.2, 0.1, -1.2, 0.4, 0.0, -0.1);
	 term(2, 5, 0, 4.2, -1.1, -4.4, -34.6, -0.2, 0.1);
	 term(2, 4, 0, 8.6, -2.9, -33.4, -97.0, 0.2, 0.1);
	 term(3, 5, 0, 0.1, -0.2, 2.1, -1.2, 0.0, 0.1);
	 term(3, 4, 0, -4.6, 9.3, 38.2, 19.8, 0.1, 0.1);
	 term(3, 3, 0, -0.5, 1.7, 23.5, 7.0, 0.0, 0.0);
	 term(4, 4, 0, 0.2, 0.8, 3.3, -1.5, -0.2, -0.1);
	 term(4, 3, 0, 0.9, 1.7, 17.9, -9.1, -0.1, 0.0);
	 term(4, 2, 0, -0.4, -0.4, -6.2, 4.8, 0.0, 0.0);
	 term(5, 3, 0, -1.6, -0.5, -2.2, 7.0, 0.0, 0.0);
	 term(5, 2, 0, -0.4, -0.1, -0.7, 5.5, 0.0, 0.0);
	 term(5, 1, 0, 0.2, 0.0, 0.0, -3.5, 0.0, 0.0);
	 term(6, 2, 0, -0.3, 0.2, 2.1, 2.7, 0.0, 0.0);
	 term(6, 1, 0, 0.1, -0.1, -1.4, -1.4, 0.0, 0.0);
	 term(6, 0, 0, -0.1, 0.1, 1.4, 0.7, 0.0, 0.0);

	 l = p2 * frac(0.1254046+m8/p2+((4982.8-21.3*t)*t+dl)/1296.0e3);
	 b = (54.77+(0.26+0.06*t)*t+db)/arc;
	 r = 30.072984+(0.001234+0.000003*t)*t+dr*1.0e-5;

	 dl = 1.04+0.02*cos(m8); dr = 0.27*sin(m8); db = 0.03*sin(m8);

	 posvel();

	return rp;
 }

/*--------------------- Pluto ------------------------------------------*/

Vec3 Plan200::Pluto (double t)   // position of Pluto at time t
 {
 /*
	 Ecliptic coordinates rp (in A.U.) and velocity vp (in A.U./day)
	 of Pluto for Equinox of Date.
	 t is the time in Julian centuries since J2000.
  */

	 const double arc = 206264.81;    // 3600*180/pi = ''/r
	 const double p2 = M_PI * 2.0;

	 int i;

	 tt = t;
	 dl = 0.0; dr = 0.0; db = 0.0;
	 m5 = p2 * frac(0.0565314+8.4302963*t);
	 m6 = p2 * frac(0.8829867+3.3947688*t);
	 m1 = p2 * frac(0.0385795+0.4026667*t);
	 c3[0] = 1.0; s3[0] = 0.0; c3[1] = cos(m1); s3[1] = sin(m1);
	 for (i=2; i<7; i++)
		addthe(c3[i-1],s3[i-1],c3[1],s3[1],c3[i],s3[i]);

	 // Kepler terms and perturbations by Jupiter
	 c[3] = 1.0; s[3] = 0.0; c[4] = cos(m5); s[4] = sin(m5);
	 for (i=3; i>1; i--)
		addthe(c[i],s[i],c[4],-s[4],c[i-1],s[i-1]);
	 addthe(c[4],s[4],c[4],s[4],c[5],s[5]);
	 term(1, 3, 0, 0.06,100924.08, -960396.0,15965.1, 51987.68,-24288.76);
	 term(2, 3, 0, 3274.74,17835.12, -118252.2,3632.4, 12687.49,-6049.72);
	 term(3, 3, 0, 1543.52,4631.99, -21446.6,1167.0, 3504.0,-1853.1);
	 term(4, 3, 0, 688.99,1227.08, -4823.4,213.5, 1048.19,-648.26);
	 term(5, 3, 0, 242.27, 415.93, -1075.4, 140.6, 302.33, -209.76);
	 term(6, 3, 0, 138.41, 110.91, -308.8, -55.3, 109.52, -93.82);
	 term(3, 2, 0, -0.99, 5.06, -25.6, 19.8, 1.26, -1.96);
	 term(2, 2, 0, 7.15, 5.61, -96.7, 57.2, 1.64, -2.16);
	 term(1, 2, 0, 10.79, 23.13, -390.4, 236.4, -0.33, 0.86);
	 term(0, 4, 0, -0.23, 4.43, 102.8, 63.2, 3.15, 0.34);
	 term(1, 4, 0, -1.1, -0.92, 11.8, -2.3, 0.43, 0.14);
	 term(2, 4, 0, 0.62, 0.84, 2.3, 0.7, 0.05, -0.04);
	 term(3, 4, 0, -0.38, -0.45, 1.2, -0.8, 0.04, 0.05);
	 term(4, 4, 0, 0.17, 0.25, 0.0, 0.2, -0.01, -0.01);
	 term(3, 1, 0, 0.06, 0.07, -0.6, 0.3, 0.03, -0.03);
	 term(2, 1, 0, 0.13, 0.2, -2.2, 1.5, 0.03, -0.07);
	 term(1, 1, 0, 0.32, 0.49, -9.4, 5.7, -0.01, 0.03);
	 term(0, 1, 0, -0.04, -0.07, 2.6, -1.5, 0.07, -0.02);

	 // perturbations by Saturn
	 c[4] = cos(m6); s[4] = sin(m6);
	 for (i=3; i>1; i--)
		addthe(c[i],s[i],c[4],-s[4],c[i-1],s[i-1]);
	 term(1, 2, 0, -29.47, 75.97, -106.4, -204.9, -40.71, -17.55);
	 term(0, 4, 0, -13.88, 18.2, 42.6, -46.1, 1.13, 0.43);
	 term(1, 4, 0, 5.81, -23.48, 15.0, -6.8, -7.48, 3.07);
	 term(2, 4, 0, -10.27, 14.16, -7.9, 0.4, 2.43, -0.09);
	 term(3, 4, 0, 6.86, -10.66, 7.3, -0.3, -2.25, 0.69);
	 term(2, 1, 0, 4.23, 2.0, 0.0, -2.2, -0.24, 0.12);
	 term(1, 1, 0, -5.04, -0.83, -9.2, -3.1, 0.79, -0.24);
	 term(0, 1, 0, 4.25, 2.48, -5.9, -3.3, 0.58, 0.02);

	 // perturbations by Jupiter and Saturn
	 double x, y;
	 y = (m5-m6); x = cos(y); y = sin(y);
	 dl = dl-9.11*x+0.12*y; dr = dr-3.4*x-3.3*y; db = db+0.81*x+0.78*y;
	 addthe(x,y,c3[1],s3[1],x,y);
	 dl = dl+5.92*x+0.25*y; dr = dr+2.3*x-3.8*y; db = db-0.67*x-0.51*y;

	 l = p2 * frac(0.6232469+m1/p2+dl/1296.0e3);
	 b = -6.8232495189e-2 + db/arc;
	 r = 40.7247248+dr*1.0e-5;

	 dl = 0.69+0.34*cos(m1)+0.12*cos(2*m1)+0.05*cos(3*m1);
	 dr = 6.66*sin(m1)+1.64*sin(2*m1);
	 db = -0.08*cos(m1)-0.17*sin(m1)-0.09*sin(2*m1);

	 // position
	 double cl, sl, cb, sb;
	 Mat3 mx;

	 cl = cos(l); sl = sin(l); cb = cos(b); sb = sin(b);
	 rp[0] = r*cl*cb; rp[1] = r*sl*cb; rp[2] = r*sb;
	 mx = pmatecl (-0.5000021, t);   // 1950.0 -> Equinox of Date
	 rp = mxvct (mx, rp);

	 // velocity
	 vp[0] = (dr*cl*cb - dl*r*sl*cb - db*r*cl*sb) * 1E-4;
	 vp[1] = (dr*sl*cb + dl*r*cl*cb - db*r*sl*sb) * 1E-4;
	 vp[2] = (dr*sb                 + db*r*cb) * 1E-4;

	return rp;
 }

/***************************************************************************/
/* =========================================================================
  Procedures for calculating the positions of various moons in the
  solar system. The algorithms for the procedures of this unit are
  taken from the Explanatory Supplement to the Astronomical Almanac
  (edited by Kenneth Seidelmann, University Science Books, Mill Valley,
  California, 1992. A number of errors in chapter 6 of this book have
  been identified and corrected in the code used here.

  Copyright :Gerhard HOLTKAMP          25-MAR-2014   
  ========================================================================= */


/*--------------------- MarPhobos -----------------------------------------*/
void MarPhobos (double t, Vec3& rs, Vec3& vs)
 {
	/*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of the Mars moon Phobos with respect to Mars for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double p2 = M_PI * 2.0;
	 const double ax =6.26974E-5;
	 const double nm = 1128.844556;
	 const double ee = 0.0150;
	 const double gam = 1.10*M_PI/180.0;

	 double thet, ll, pp, na, ja;
	 double al, dy, yr;
	 Mat3 mt1;

	 // time from 11-NOV-1971
	 dy = t * 36525.0 + 10278.5;
	 yr = dy / 365.25;

	 // normalized angles
	 thet = p2 * frac ((327.9 - 0.43533*dy) / 360.0);
	 ll = p2 * frac ((232.41 + nm*dy + 0.00124*yr*yr) / 360.0);
	 pp = p2 * frac ((278.96 + 0.43526*dy) / 360.0);
	 na = p2 * frac ((47.39 - 0.0014*yr) / 360.0);
	 ja = p2 * frac ((37.27 + 0.0008*yr) / 360.0);

	 // in-orbit-plane position vector
	 al = ll - pp; // mean anomaly
	 al = eccanom (al, ee); // eccentric anomaly
	 rs[0] = ax * (cos(al) - ee);
	 rs[1] = ax * sqrt(1.0 - ee*ee) * sin(al);
	 rs[2] = 0.0;

	 // in-orbit-plane velocity vector
	 dy = cos(al);
	 yr = 1.0 - ee*dy;
	 yr = 1.235083014E-3 / yr;
	 vs.assign(-yr*sin(al),yr*sqrt(1.0-ee*ee)*dy,0.0);

	 // convert to equatorial coordinates
	 al = pp - thet - na;
	 mt1 = zrot (-al);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = xrot (-gam);
	 rs =	mxvct (mt1, rs);
	 vs = mxvct (mt1, vs);
	 mt1 = zrot (-thet);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = xrot (-ja);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = zrot (-na);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

	 // convert to ecliptic coordinates of date
	 mt1 = pmatequ (-0.500002096, t);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 rs =	equecl (t, rs);
	 vs =	equecl (t, vs);

 }

/*--------------------- MarDeimos ------------------------------------------*/
void MarDeimos (double t, Vec3& rs, Vec3& vs)
 {
	/*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of the Mars moon Deimos with respect to Mars for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================== */

	 const double p2 = M_PI * 2.0;
	 const double ax =1.56828E-4;
	 const double nm = 285.161888;
	 const double ee = 0.0004;
	 const double gam = 1.79*M_PI/180.0;

	 double thet, ll, pp, na, ja;
	 double al, dy, yr;
	 Mat3 mt1;

	 // time from 11-NOV-1971
	 dy = t * 36525.0 + 10278.5;
	 yr = dy / 365.25;

	 // normalized angles
	 thet = p2 * frac ((240.38 - 0.01801*dy) / 360.0);
	 ll = p2 * frac ((196.55 - 0.01801*dy) / 360.0);
	 ll = p2 * frac ((28.96 + nm*dy - 0.27*sin(ll)) / 360.0);
	 pp = p2 * frac ((111.7 + 0.01798*dy) / 360.0);
	 na = p2 * frac ((46.37 - 0.0014*yr) / 360.0);
	 ja = p2 * frac ((36.62 + 0.0008*yr) / 360.0);

	 // in-orbit- plane position vector
	 al = ll - pp; // mean anomaly
	 al = eccanom (al, ee); // eccentric anomaly
	 rs[0] = ax * (cos(al) - ee);
	 rs[1] = ax * sqrt(1.0 - ee*ee) * sin(al);
	 rs[2] = 0.0;

	 // in-orbit-plane velocity vector
	 dy = cos(al);
	 yr = 1.0 - ee*dy;
	 yr = 7.8046400669E-4 / yr;
	 vs.assign(-yr*sin(al),yr*sqrt(1.0-ee*ee)*dy,0.0);

	 // convert to equatorial coordinates
	 al = pp - thet - na;
	 mt1 = zrot (-al);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = xrot (-gam);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = zrot (-thet);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = xrot (-ja);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = zrot (-na);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

	 // convert to ecliptic coordinates of date
	 mt1 = pmatequ (-0.500002096, t);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 rs =	equecl (t, rs);
	 vs =	equecl (t, vs);

 }

/*--------------------- PosJIo ------------------------------------------*/
Vec3 PosJIo (double t)
 {
  /* ---------------------------------------------------------------------- }
	 Ecliptic coordinates (in A.U.) rs
	 of Io with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ======================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ax = 0.002819347;
	 const double omj = 99.99754*M_PI/180.0;
	 const double jj = 1.30691*M_PI/180.0;
	 const double ii = 3.10401*M_PI/180.0;

	 double l1, l2, phl, p1, p2, p3, p4;
	 double pij, th1, th2, psi, gg;
	 double xi, nu, ze;
	 double d;  // : EXTENDED;
	 Vec3 rs;
	 Mat3 mt1;

	 // general jupiter moon data
	 d = t*36525.0 + 8544.5;
	 l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
	 l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
	 phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
	 p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
	 p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
	 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
	 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
	 pij = 0.23510274429; // 13.470395
	 th1 = pi2* frac((308.365749 - 0.1328061*d)/360.0);
	 th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
	 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
	 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);

	 // io specific lines
	 xi = -41279.0 * cos(2.0*(l1-l2))*1E-7;
	 nu = -5596.0*sin(p3-p4) - 2198.0*sin(p1+p2-2.0*(pij+gg))
				+1321.0*sin(phl) - 1157.0*sin(l1-2.0*l2+p4)
			-1940.0*sin(l1-2.0*l2+p3) + 791.0*sin(l1-2.0*l2+p2)
				+ 82363.0*sin(2.0*(l1-l2));
	 nu *= 1E-7;
	 ze = 7038.0*sin(l1-th1+nu) + 1835.0*sin(l1-th2+nu);
	 ze *= 1E-7;

	 // convert to rectangular coordinates in Jovian equatorial

	 rs.assign (ax * (1+xi) * cos(l1-psi+nu),
					ax * (1+xi) * sin(l1-psi+nu),
					ax * ze);

	 // convert into ecliptic of 1950.0
	 mt1 = xrot (-ii);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-(psi-omj));
	 rs =	mxvct (mt1, rs);
	 mt1 = xrot (-jj);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-omj);
	 rs =	mxvct (mt1, rs);

	 // convert into ecliptic of date
	 mt1 = pmatecl (-0.500002096, t);
	 rs =	mxvct (mt1, rs);

	return rs;

 }

/*--------------------- PosEuropa ------------------------------------------*/
Vec3 PosEuropa (double t)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs
	 of Europa with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 =====================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ax = 0.004485872;
	 const double omj = 99.99754*M_PI/180.0;
	 const double jj = 1.30691*M_PI/180.0;
	 const double ii = 3.10401*M_PI/180.0;

	 double l1, l2, l3, phl, p2, p3, p4;
	 double pij, th2, th3, psi, gg;
	 double xi, nu, ze;
	 double d; // : EXTENDED;
	 Vec3 rs;
	 Mat3 mt1;


	 // general jupiter moon data
	 d = t*36525.0 + 8544.5;
	 l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
	 l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
	 l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
	 phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
	 // p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
	 p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
	 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
	 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
	 pij = 0.23510274429; // 13.470395
	 th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
	 th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
	 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
	 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);

	 // Europa specific lines
	 xi = -3187*cos(l2-p3) - 1738*cos(l2-p4)
				+93748*cos(l1-l2);
	 xi *= 1E-7;
	 nu = -1158*sin(2*(psi-pij)) + 1715*sin(-2*(pij+gg)+th3+psi)
				-1846*sin(gg) + 2397*sin(p3-p4) - 3172*sin(phl)
				-1993*sin(l2-l3) + 1844*sin(l2-p2)
				+6394*sin(l2-p3) + 3451*sin(l2-p4)
				+4159*sin(l1-2*l2+p4) + 7571*sin(l1-2*l2+p3)
				-1491*sin(l1-2*l2+p3) - 185640*sin(l1-l2)
				-803*sin(l1-l3) + 915*sin(2*(l1-l2));
	 nu *= 1E-7;
	 ze = 81575*sin(l2-th2+nu) + 4512*sin(l2-th3+nu)
				-3286*sin(l2-psi+nu);
	 ze *= 1E-7;

	 // convert to rectangular coordinates in Jovian equatorial

	 rs.assign (ax * (1+xi) * cos(l2-psi+nu),
					ax * (1+xi) * sin(l2-psi+nu),
					ax * ze);

	 // convert into ecliptic of 1950.0

	 mt1 = xrot (-ii);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-(psi-omj));
	 rs =	mxvct (mt1, rs);
	 mt1 = xrot (-jj);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-omj);
	 rs =	mxvct (mt1, rs);

	 // convert into ecliptic of date
	 mt1 = pmatecl (-0.500002096, t);
	 rs =	mxvct (mt1, rs);

	return rs;

  }

/*--------------------- PosGanymede --------------------------------------*/
Vec3 PosGanymede (double t)
 {
  /*
	 Ecliptic coordinates (in A.U.)
	 of Ganymede with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 =====================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ax = 0.007155352;
	 const double omj = 99.99754*M_PI/180.0;
	 const double jj = 1.30691*M_PI/180.0;
	 const double ii = 3.10401*M_PI/180.0;

	 double l1, l2, l3, l4, phl, p1, p2, p3, p4;
	 double pij, th2, th3, th4, psi, gg;
	 double xi, nu, ze;
	 double d;  // : EXTENDED;
	 Vec3 rs;
	 Mat3 mt1;

	 // general jupiter moon data
	 d = t*36525.0 + 8544.5;
	 l1 = pi2 * frac((106.07859 + 203.4889553630643*d)/360.0);
	 l2 = pi2 * frac((175.733787 + 101.3747245566245*d)/360.0);
	 l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
	 l4 = pi2 * frac((84.455823 + 21.57107087517961*d)/360.0);
	 phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
	 p1 = pi2 * frac((82.380231 + 0.16102275*d)/360.0);
	 p2 = pi2 * frac((128.960393 + 0.04645644*d)/360.0);
	 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
	 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
	 pij = 0.23510274429; // 13.470395
	 th2 = pi2* frac((100.438938 - 0.03261535*d)/360.0);
	 th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
	 th4 = pi2* frac((322.746564 - 0.00176018*d)/360.0);
	 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
	 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);

	 // Ganymede specific lines
	 xi = -14691*cos(l3-p3) - 1758*cos(2*(l3-l4))
				+6333*cos(l2-l3);
	 xi *= 1E-7;
	 nu = -1488*sin(2*(pij+psi)) + 411*sin(-th3+psi)
				+346*sin(-th4+psi) - 2338*sin(gg)
				+6558*sin(p3-p4) + 523*sin(p1+p3-2*(pij+gg))
				+341*sin(phl) - 943*sin(l3-l4)
				+29387*sin(l3-p3) + 15800*sin(l3-p4)
				+3218*sin(2*(l3-l4)) + 226*sin(3*l3-2*l4)
				-12038*sin(l2-l3) - 662*sin(l1-2*l2+p4)
				-1246*sin(l1-2*l2+p3) + 699*sin(l1-2*l2+p2)
				+217*sin(l1-l3);
	 nu *= 1E-7;
	 ze = -2793*sin(l3-th2+nu) + 32387*sin(l3-th3+nu)
				+6871*sin(l3-th4+nu) - 16876*sin(l3-psi+nu);
	 ze *= 1E-7;

	 // convert to rectangular coordinates in Jovian equatorial

	 rs.assign (ax * (1+xi) * cos(l3-psi+nu),
					ax * (1+xi) * sin(l3-psi+nu),
					ax * ze);

	 // convert into ecliptic of 1950.0

	 mt1 = xrot (-ii);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-(psi-omj));
	 rs =	mxvct (mt1, rs);
	 mt1 = xrot (-jj);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-omj);
	 rs =	mxvct (mt1, rs);

	 // convert into ecliptic of date
	 mt1 = pmatecl (-0.500002096, t);
	 rs =	mxvct (mt1, rs);

	return rs;

 }

/*--------------------- PosCallisto --------------------------------------*/
Vec3 PosCallisto (double t)
 {
  /*
	 Ecliptic coordinates (in A.U.)
	 of Callisto with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 =====================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ax = 0.012585436;
	 const double omj = 99.99754*M_PI/180.0;
	 const double jj = 1.30691*M_PI/180.0;
	 const double ii = 3.10401*M_PI/180.0;

	 double l3, l4, p3, p4;
	 double pij, th3, th4, psi, gp, gg, ph2;
	 double xi, nu, ze;
	 double d;  // : EXTENDED;
	 Vec3 rs;
	 Mat3 mt1;

	 // general jupiter moon data
	 d = t*36525.0 + 8544.5;
	 l3 = pi2 * frac((120.5613855 + 50.31760915340462*d)/360.0);
	 l4 = pi2 * frac((84.455823 + 21.57107087517961*d)/360.0);
	 // phl = pi2* frac((184.415351 + 0.17356902*d)/360.0);
	 p3 = pi2 * frac((187.550171 + 0.00712408*d)/360.0);
	 p4 = pi2 * frac((335.309254 + 0.00183939*d)/360.0);
	 pij = 0.23510274429; // 13.470395
	 th3 = pi2* frac((118.908928 - 0.00717678*d)/360.0);
	 th4 = pi2* frac((322.746564 - 0.00176018*d)/360.0);
	 psi = pi2* frac((316.500101 - 0.00000248*d)/360.0);
	 gp = pi2 * frac((31.9785280244 + 0.033459733896*d)/360.0);
	 gg = pi2 * frac((30.2380210168 + 0.08309256178969453*d)/360.0);
	 ph2 = 0.91009489942; // 52.1445966929

	 // Callisto specific lines
	 xi = 1656*cos(l4-p3) - 73328*cos(l4-p4)
				+182*cos(l4-pij) - 541*cos(l4+p4-2*(pij+gg))
				-269*cos(2*(l4-p4)) + 974*cos(l3-l4);
	 xi *= 1E-7;
	 nu = -407*sin(2*(psi-p4)) + 309*sin(-2*p4+th4+psi)
				-4840*sin(2*(psi-pij)) + 2074*sin(-th4+psi)
				-5605*sin(gg) - 204*sin(2*gg)
				-495*sin(5*gp-2*gg+ph2) + 234*sin(p4-pij)
				-6112*sin(p3-p4) - 3318*sin(l4-p3)
				+145573*sin(l4-p4) + 178*sin(l4-pij-gg)
				-363*sin(l4-pij) + 1085*sin(l4+p4-2*(pij+gg))
				+672*sin(2*(l4-p4)) + 218*sin(2*(l4-pij-gg))
				+167*sin(2*l4-th4-psi) - 142*sin(2*(l4-psi))
				+148*sin(l3-2*l4+p4) - 390*sin(l3-l4)
				-195*sin(2*(l3-l4)) + 185*sin(3*l3-7*l4+4*p4);
	 nu *= 1E-7;
	 ze = 773*sin(l4-2*(pij+gg)+psi+nu)
				-5075*sin(l4-th3+nu) + 44300*sin(l4-th4+nu)
				-76493*sin(l4-psi+nu);
	 ze *= 1E-7;

	 rs.assign(ax * (1+xi) * cos(l4-psi+nu),
				  ax * (1+xi) * sin(l4-psi+nu),
				  ax * ze);

	 // convert into ecliptic of 1950.0

	 mt1 = xrot (-ii);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-(psi-omj));
	 rs =	mxvct (mt1, rs);
	 mt1 = xrot (-jj);
	 rs =	mxvct (mt1, rs);
	 mt1 = zrot (-omj);
	 rs =	mxvct (mt1, rs);

	 // convert into ecliptic of date
	 mt1 = pmatecl (-0.500002096, t);
	 rs =	mxvct (mt1, rs);

	return rs;

 }

/*-------------------- Mimas, Enceladus, Dione ---------------------*/
Vec3 PosSMimas (double t)
{
 /* Ecliptic coordinates (in A.U.) of Mimas with respect to Saturn 
    for Equinox of Date.
    t is the time in Julian centuries since J2000.
  ===================================================================*/
 const double ie = 28.0771*M_PI/180.0;
 const double oe = 168.8263*M_PI/180.0;

 double d, dt, tt, a, n, ecc, gam, th1, l1, p, m;
 Vec3 rs;
 Mat3 mt1; 

 d = 36525.0*t + 40452.0;
 dt = d / 365.25;
 tt = 5.0616*(((36525*t + 18262.577000)/365.25)+1950.0-1866.06);
 tt *= M_PI / 180.0;
 a = 0.00124171;
 n = 381.994516;
 ecc = 0.01986;
 gam = 0.0274017; // 1.570;
 th1 = 49.4 - 365.025*dt;
 l1 = 128.839 + n*d - 43.415*sin(tt) - 0.714*sin(3.0*tt) - 0.020*sin(5.0*tt);
 p = 107.0 + 365.560*dt;

 // in-orbit-plane position vector
 m = l1 - p; // mean anomaly
 m = m * M_PI / 180.0;
 m = eccanom (m, ecc);  // eccentric anomaly
 rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);

 // convert to ecliptic coordinates
 m = p - th1;
 m = m * M_PI / 180.0;
 mt1 = zrot (-m);
 rs = mxvct (mt1, rs);
 mt1 = xrot (-gam);
 rs = mxvct (mt1, rs);
 m = th1*M_PI/180.0 - oe;
 mt1 = zrot (-m);
 rs = mxvct (mt1, rs);
 mt1 = xrot (-ie);
 rs = mxvct (mt1, rs);
 mt1 = zrot (-oe);
 rs = mxvct (mt1, rs);

 // convert to ecliptic coordinates of date
 mt1 = pmatecl (-0.500002096, t);
 rs = mxvct (mt1, rs);

 return rs;
}

Vec3 PosSEnceladus (double t)
{
 /* Ecliptic coordinates (in A.U.) of Enceladus with respect to Saturn 
    for Equinox of Date.
    t is the time in Julian centuries since J2000.
  ===================================================================*/
 const double ie = 28.0771*M_PI/180.0;
 const double oe = 168.8263*M_PI/180.0;

 double d, dt, /*tt,*/ a, n, ecc, gam, th1, l1, p, m;
 Vec3 rs;
 Mat3 mt1;

 d = 36525.0*t + 40452.0;
 dt = d / 365.25;
 a = 0.00158935;
 n = 262.7319052;
 ecc = 0.00532;
 gam = 0.000628319; //  0.036;
 th1 = 145.0 - 152.7*dt;
 l1 = 200.155 + n*d + 0.256333 * sin((59.4+32.73*dt)*M_PI/180.0)
      + 0.217333 * sin((119.2+93.18*dt)*M_PI/180.0); 
 p = 312.7 + 123.42*dt;

 // in-orbit-plane position vector
 m = l1 - p; // mean anomaly
 m = m * M_PI / 180.0;
 m = eccanom (m, ecc);  // eccentric anomaly
 rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);

 // convert to ecliptic coordinates
 m = p - th1;
 m = m * M_PI / 180.0;
 mt1 = zrot (-m);
 rs = mxvct (mt1, rs);
 mt1 = xrot (-gam);
 rs = mxvct (mt1, rs);
 m = th1*M_PI/180.0 - oe;
 mt1 = zrot (-m);
 rs = mxvct (mt1, rs);
 mt1 = xrot (-ie);
 rs = mxvct (mt1, rs);
 mt1 = zrot (-oe);
 rs = mxvct (mt1, rs);

 // convert to ecliptic coordinates of date
 mt1 = pmatecl (-0.500002096, t);
 rs = mxvct (mt1, rs);

 return rs;
}
 
Vec3 PosSDione (double t)
{
 /* Ecliptic coordinates (in A.U.) of Dione with respect to Saturn 
    for Equinox of Date.
    t is the time in Julian centuries since J2000.
  ===================================================================*/

 const double ie = 28.0771*M_PI/180.0;
 const double oe = 168.8263*M_PI/180.0;

 double d, dt, /*tt,*/ a, n, ecc, gam, th1, l1, p, m;
 Vec3 rs;
 Mat3 mt1;

 d = 36525.0*t + 40452.0;
 dt = d / 365.25;
 a = 0.00252413;
 n = 131.534920026;
 ecc = 0.001715;
 gam = 0.0005044; // 0.0289;
 th1 = 228.0 - 30.6197 * dt;
 l1 = 255.1183 + n*d - 0.0146667 * sin((59.4+32.73*dt)*M_PI/180.0)
      - 0.0125 * sin((119.2+93.18*dt)*M_PI/180.0);
 p = 173.6 + 30.8381 * dt;

 // in-orbit-plane position vector
 m = l1 - p; // mean anomaly
 m = m * M_PI / 180.0;
 m = eccanom (m, ecc);  // eccentric anomaly
 rs.assign (a*(cos(m) - ecc), a*sqrt(1.0 - ecc*ecc)*sin(m), 0);

 // convert to ecliptic coordinates
 m = p - th1;
 m = m * M_PI / 180.0;
 mt1 = zrot (-m);
 rs = mxvct (mt1, rs);
 mt1 = xrot (-gam);
 rs = mxvct (mt1, rs);
 m = th1*M_PI/180.0 - oe;
 mt1 = zrot (-m);
 rs = mxvct (mt1, rs);
 mt1 = xrot (-ie);
 rs = mxvct (mt1, rs);
 mt1 = zrot (-oe);
 rs = mxvct (mt1, rs);

 // convert to ecliptic coordinates of date
 mt1 = pmatecl (-0.500002096, t);
 rs = mxvct (mt1, rs);

 return rs;
}

/*--------------------- JupIo --------------------------------------*/
void JupIo (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Io with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double dlt = 3.472222E-4;  // 30 sec

	 double td;

	 td = dlt / 36525;
	 rs = PosJIo (t-td);
	 vs = PosJIo (t+td);
	 vs -= rs;
	 td = 0.5/dlt;
	 vs *= td;
	 rs = PosJIo (t);

  }

/*--------------------- JupEuropa --------------------------------------*/
void JupEuropa (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Europa with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double dlt = 3.472222E-4;  // 30 sec

	 double td;

	 td = dlt / 36525;
	 rs = PosEuropa (t-td);
	 vs = PosEuropa (t+td);
	 vs -= rs;
	 td = 0.5/dlt;
	 vs *= td;
	 rs = PosEuropa (t);

  }

/*--------------------- JupGanymede ------------------------------------*/
void JupGanymede (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Ganimede with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double dlt = 3.472222E-4;  // 30 sec

	 double td;

	 td = dlt / 36525;
	 rs = PosGanymede (t-td);
	 vs = PosGanymede (t+td);
	 vs -= rs;
	 td = 0.5/dlt;
	 vs *= td;
	 rs = PosGanymede (t);

  }

/*--------------------- JupCallisto ------------------------------------*/
void JupCallisto (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Callisto with respect to Jupiter for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double dlt = 3.472222E-4;  // 30 sec

	 double td;

	 td = dlt / 36525;
	 rs = PosCallisto (t-td);
	 vs = PosCallisto (t+td);
	 vs -= rs;
	 td = 0.5/dlt;
	 vs *= td;
	 rs = PosCallisto (t);

  }

/*----------------------- SatRhea ------------------------------------*/
void SatRhea (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Rhea with respect to Saturn for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ie = 28.0771*M_PI/180.0;
	 const double oe = 168.8263*M_PI/180.0;
	 const double kap = 57.29578*M_PI/180.0;
	 const double ax = 0.003524;
	 const double sg0 = 5.7682811893E-3;  // sin(0.3305)

	 double d, td, tt, pp, ot, nt;
	 double lam, ii, omg, omb, ee;
	 Mat3 mt1;

	 // various times
	 tt = t*36525.0;
	 d = tt + 40452.0;
	 td = 0.5219*(tt+40177)/365.25;
	 tt = d / 365.25;

	 // Titan perturbations
	 pp = pi2 * frac((305.0 + 10.2077*tt)/360.0);
	 ot = pi2 * frac((276.49 + td)/360.0);
	 nt = pi2 * frac((44.5 - td)/360.0);

	 // osculating orbit elements
	 ee = 0.00021*sin(pp) + 0.001*sin(ot);
	 omg = 0.00021*cos(pp) + 0.001*cos(ot);
	 omb = atan2(ee,omg);
	 ii = sin(omb);
	 if (fabs(ii) > 0.5) ee = fabs(ee/ii);
	 else ee = fabs(omg/cos(omb));
	 pp = pi2 * frac((356.87 - 10.2077*tt)/360.0);
	 ot = sin(pp);
	 lam = pi2 * frac((359.4727 + 79.69004007*d)/360.0);
	 lam = lam + kap*sg0*tan(0.5*ie)*ot;
	 ii = ie - 7.9412480966E-4 + kap*sg0*cos(pp)
				+ 3.5081117965E-4*cos(nt);
	 omg = oe - 1.3613568166E-4
				 + (kap*sg0*ot + 3.5081117965E-4*sin(nt))/sin(ie);

	 // in-orbit-plane position vector
	 pp = lam - omb; // mean anomaly
	 pp = eccanom (pp, ee);  // eccentric anomaly
	 rs.assign (ax * (cos(pp) - ee), ax * sqrt(1.0 - ee*ee) * sin(pp), 0);

	 // in-orbit-plane velocity vector
	 d = cos(pp);
	 td = 1.0 - ee*d;
	 td = 4.9000413575E-3 / td;
	 vs.assign (-td*sin(pp), td*sqrt(1.0-ee*ee)*d, 0);

	 // convert to ecliptic coordinates
	 pp = omb - omg;
	 mt1 = zrot (-pp);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = xrot (-ii);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = zrot (-omg);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

	 // convert to ecliptic coordinates of date
	 mt1 = pmatecl (-0.500002096, t);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

  }

/*----------------------- SatTitan ------------------------------------*/
void SatTitan (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Titan with respect to Saturn for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ===================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ie = 28.0771*M_PI/180.0;
	 const double oe = 168.8263*M_PI/180.0;
	 const double kap = 57.29578*M_PI/180.0;
	 const double ax = 0.00816765;
	 const double sg0 = 5.7682811893E-3;  // sin(0.3305)

	 double d, tt, ts, gg, ls, is, os, lms, psi;
	 double lam, ii, omg, omb, ee, s4;
	 Mat3 mt1;

	 // various times
	 ts = t + 1.0;
	 tt = t*36525.0;
	 d = tt+40177;
	 tt = d / 365.25;

	 // solar perturbations
	 ls = pi2 * frac((175.4762 + 1221.5515*ts)/360.0);
	 is = pi2 * frac((2.4891 + 0.002435*ts)/360.0);
	 os = pi2 * frac((113.35 - 0.2597*ts)/360.0);
	 lms = pi2* frac((267.2635 + 1222.1136*ts)/360.0);

	 gg = pi2 * frac((41.28 - 0.5219*tt)/360.0);
	 ii = ie - 1.0828022679E-2 + kap*5.2185107774E-3*cos(gg);
	 s4 = sin(gg);
	 omg = oe - 2.4748768793E-3 + kap*5.2185107774E-3*s4/sin(ie);
	 omb = pi2 * frac((275.837 + 0.5219*tt)/360.0);

	 psi = sin(is)*sin(omg-os);
	 gg = cos(omg-os);
	 ts = cos(is)*sin(ii) - sin(is)*cos(ii)*gg;
	 psi = atan2(psi,ts);
	 ts = cos(is)*sin(ii)*gg - sin(is)*cos(ii);
	 gg = sin(ii)*sin(omg-os);
	 gg = atan2(gg,ts);
	 lms = lms - gg - os;
	 gg = omb - omg - psi;

	 // osculating orbit elements
	 ts = 2*(lms-gg);
	 ee = 0.028815 - 0.000184*cos(2*gg) + 0.000073*cos(ts);
	 omb = omb + kap*(0.00630*sin(2*gg) + 0.0025*sin(ts));
	 gg = 2*lms + psi;
	 lam = pi2 * frac((261.3121 + 22.57697385*d)/360.0)
				 +kap*(sg0*tan(0.5*ie)*s4 - 0.000176*sin(ls)
				 -0.000215*sin(2*lms) + 0.000057*sin(gg));
	 ii = ii + 0.000232*kap*cos(gg);
	 omg = omg + 0.000503*kap*sin(gg);

	 // in-orbit-plane position vector
	 gg = lam - omb; // mean anomaly
	 if (gg < 0) gg += pi2;
	 gg = eccanom (gg, ee); // eccentric anomaly
	 rs.assign (ax * (cos(gg) - ee), ax * sqrt(1.0 - ee*ee) * sin(gg), 0);

	 // in-orbit-plane velocity vector
	 d = cos(gg);
	 ts = 1.0 - ee*d;
	 ts = 3.2183196288E-3 / ts;
	 vs.assign (-ts*sin(gg), ts*sqrt(1.0-ee*ee)*d, 0);

	 // convert to ecliptic coordinates
	 gg = omb - omg;
	 mt1 = zrot (-gg);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = xrot (-ii);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = zrot (-omg);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

	 // convert to ecliptic coordinates of date
	 mt1 = pmatecl (-0.500002096, t);
	 rs = mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

  }


/*----------------------- NepTriton ------------------------------------*/
void NepTriton (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Triton with respect to Neptune for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ax = 0.00237142;
	 const double gam = 158.996*M_PI/180.0;

	 double d, nn, u, jj, ap;
	 Mat3 mt1;

	 // days since epoch
	 d = t*36525 + 18262.5;

	 // pole of fixed plane in 1950.0
	 nn = pi2 * frac((359.28 + 54.308*t)/360.0);
	 rs.assign (1.0,
					pi2 * frac((298.72+2.58*sin(nn)-0.04*sin(2*nn))/360.0),
					pi2 * frac((42.63-1.9*cos(nn)+0.01*cos(2*nn))/360.0));
	 rs =	polcar (rs);
	 mt1 = pmatequ (0.0, -0.500002096);
	 rs =	mxvct (mt1, rs);
	 rs =	carpol (rs);
	 jj = M_PI/2.0;
	 ap = rs[1] + jj;   // R.A. of ascending node
	 jj = jj - rs[2];   // inclination

	 // in-orbit-plane position vector

	 u = pi2 * frac((200.913 + 61.2588532*d)/360.0);
	 rs.assign (ax * cos(u), ax * sin(u), 0);

	 // velocity is perpendicular to radius vector for e=0
	 u = 1.0 / abs(rs);
	 vs = u * rs;
	 mt1 = zrot (-M_PI/2.0);
	 vs =	mxvct (mt1, vs);
	 u = 2.53538612E-3;  // mean orbital speed
	 vs *= u;

	 // convert to equatorial coordinates
	 mt1 = xrot (-gam);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 nn = pi2 * frac((151.401 + 0.57806*d/365.25)/360.0);
	 mt1 = zrot (-nn);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = xrot (-jj);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = zrot (-ap);
	 rs = mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

	 // convert to ecliptic coordinates of date
	 mt1 = pmatequ (-0.500002096, t);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 rs =	equecl (t, rs);
	 vs =	equecl (t, vs);

  }


/*----------------------- PluCharon ------------------------------------*/
void PluCharon (double t, Vec3& rs, Vec3& vs)
 {
  /*
	 Ecliptic coordinates (in A.U.) rs and velocity vs (in A.U./day)
	 of Charon with respect to Pluto for Equinox of Date.
	 t is the time in Julian centuries since J2000.
	 ====================================================================*/

	 const double pi2 = M_PI * 2.0;
	 const double ax = 0.00013102;
	 const double jj = 94.3*M_PI/180.0;
	 const double nn = 223.7*M_PI/180.0;

	 double d, u;
	 Mat3 mt1;

	 // days since epoch
	 d = t*36525 + 6544.5;

	 // in-orbit-plane position vector
	 u = pi2 * frac((78.6 + 56.3625*d)/360.0);
	 rs.assign (ax * cos(u), ax * sin(u), 0);

	 // velocity is perpendicular to radius vector for e=0
	 d = 1.0 / abs(rs);
	 vs = d * rs;
	 mt1 = zrot (-M_PI/2.0);
	 vs =  mxvct (mt1, vs);
	 d = 1.288837578E-4;  // mean orbital speed
	 vs *= d;

	 // convert to equatorial coordinates
	 mt1 = xrot (-jj);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 mt1 = zrot (-nn);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);

	 // convert to ecliptic coordinates of date
	 mt1 = pmatequ (-0.500002096, t);
	 rs =	mxvct (mt1, rs);
	 vs =	mxvct (mt1, vs);
	 rs =	equecl (t, rs);
	 vs =	equecl (t, vs);

  }

